      program frisco
!	dust + gas radiative transfer in a circumstellar disk
!	with nested r-theta grids tranformed onto rectangular grids
!	for the rad xfer
!	expand the innner grid by 8 to resolve the inner rim better
!	start with Chiang & Goldreich disk model
!	with hydro equil for vertical structure
!	then do full dust rad xfer
!	this version calculates rotational lines with populations calculated
!	using a critical density approach, ignoring radiative excitation
!	probably could incorporate radiation into critical density
!	maybe with just one iteration

        implicit real*8 (a-h,o-z)
#include "mpif.h"
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        parameter (mw = 1024)
        parameter (mc = 256)
        parameter (mrz = 512)
        parameter (mmap = 1024)
        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /heat/ gmmdot,acom,pxr,axr,bxr,epspe,agas,sigma, &
                colin,colup,xco,xwat
        common /angls/ dcost,dcos9,rot,irot
        common /parmz/ rgrid1,r0,r1,sigma0,sigma1,ocol,pcol,qcol, &
                avisc,xmstar,xlstar,tuv,xuv,tsubl,gsubl,tburn,xburn, &
                coswind,xc,xsi,ac,asi,enhanc,xmdot,wmdot,vwind, &
                xco0,xwat0,tfrco,tfrwat,xmol,tfrmol,tdiss,coldiss, &
                kit,lit,mit,nnt,lx,lrot,reuse,lte
        common /xj/ xju,xjp,xjq,xjr,xk
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg), &
                ff(mu,mt,mr,mg),fn(mt,mr,mg)
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu), &
                phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg), &
                opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg), &
                rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml), &
                wl(ml),dl(ml,mt,mr,mg),sl(ml,mt,mr,mg),opl(ml,mt,mr,mg)
        common /vels/ wr(mt,mr,mg),wth(mt,mr,mg),wph(mt,mr,mg)
        common /taus/ taucont,tauline,taucent,iuv,iusil,iufl,iwterm
        common /cube/ cube(mc,mc,mc),ngc
        common /exptab/ expi(10239),too
        dimension mpistatus(MPI_STATUS_SIZE)
        dimension taud(mu,mt),dtac4(mt)
        dimension colr(mt,mr,mg),colz(mt,mr,mg),tauv(mt,mr,mg)
        dimension dt(mt,mr,mg),dm(mt,mr,mg)
        dimension tgz(mrz,mrz,mg),dgz(mrz,mrz,mg),ddz(mrz,mrz,mg)
        dimension fw(mw,ml,2*mt),fwj(mw)
        dimension spec(mw,mu)
        dimension frame(mc,mc),picture(mt,mr,mg)
        dimension con(mmap)
        logical ifmpi,doplot,doprint,wait
        logical bipolar,smooth,lte,reuse
        logical ifrim,inwind
        logical badnum
        character*16 monfile,pfile,qfile,rfile,hfile,ifile,lockline
        character*8 yn

!	physical constants

        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (three = 3.0d+00)
        parameter (four = 4.0+00)
        parameter (ten = 1.0d+01)
        parameter (hun = 1.0d+02)
        parameter (thrq = 0.75d+00)
        parameter (half = 0.5d+00)
        parameter (quart = 0.25d+00)
        parameter (fifth = 0.2d+00)
        parameter (tiny = 1.0d-12)
        parameter (zilch = 1.0d-200)
        parameter (pi = 3.14159265d+00)
        parameter (clt = 2.9979d+10)
        parameter (hpl = 6.6262d-27)
        parameter (hc2 = 5.955d-06)
        parameter (bk = 1.3807d-16)
        parameter (sigsb = 5.670d-05)
        parameter (au = 1.496d+13)
        parameter (rsun = 6.96d+10)
        parameter (xlsun = 3.90d+33)
        parameter (xmsun = 1.99d+33)
        parameter (secyr = 3.156d+07)
!	mass / gas particle
        parameter (xmg = 4.0d-24)
        parameter (vkms = 1.0d+05)

        sqrt2 = sqrt(two)
        alog2 = log(two)
        sqrtpi = sqrt(pi)
        tpi = two*pi
        fpi = four*pi
        third = one/3.0d+00
        hc = hpl*clt
        hck = hc/bk

        ifmpi = .true.
        if(ifmpi) then
          call mpi_init(ierr)
          call mpi_comm_size(MPI_COMM_WORLD,mysize,ierr)
          call mpi_comm_rank(MPI_COMM_WORLD,myrank,ierr)
          ifmpi = (mysize.gt.1)
        else
          mysize = 1
          myrank = 0
        endif
        
        doprint = (myrank.eq.0)
        doplot = (myrank.eq.(mysize/2))
        if(doprint) then
          print '("Running frisco with",i4," processors")', mysize
        endif
        wait = .true.
!	reuse old dust disk model?
        reuse = .false.

!	initialize expi arrays
        x = expinit()

!	initialize plot files

        if(doprint) then
          print '(" Enter name of output file ",$)'
          read '(a16)', monfile
          if(index(monfile,' ').lt.2) monfile = 'specs'
          last = index(monfile,' ')-1
          last = min(last,14)
          ifile = monfile(1:last)//'.log'
          open(unit=16,file=ifile,access='sequential')
          open(unit=12,file='data.spec',access='sequential')
          do iu = 1,64
            read (12,*), ju,wavel,flux
            spec(iu,9) = slog10(flux)
            if(ju.ne.iu) stop
          enddo
          close(unit=12)
          open(unit=12,file='data.star',access='sequential')
          do iu = 1,64
            read (12,*), ju,wavel,spec(iu,15)
            if(ju.ne.iu) stop
          enddo
          close(unit=12)
        endif
        if(ifmpi) then
          call mpi_bcast(monfile,16,MPI_CHARACTER, &
            0,MPI_COMM_WORLD,ierr)
          call mpi_bcast(spec(1,9),64,MPI_DOUBLE_PRECISION, &
            0,MPI_COMM_WORLD,ierr)
          call mpi_bcast(spec(1,15),64,MPI_DOUBLE_PRECISION, &
            0,MPI_COMM_WORLD,ierr)
        endif
        if(doplot) then
          open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
            access='sequential')
          write(8,'("empty")')
          close(unit=8)
          last = index(monfile,' ')-1
          last = min(last,14)
          rfile = monfile(1:last)//'.rot'
          open(unit=11,file=rfile,access='sequential')
          monfile = ''
        endif

!	molecular constants
!	some left over from pisco

!	J=1-0 A coefficients will be converted to band strengths
!	to use vibrational formulas
!	grot is collision rate coeff / T^1/2
!	SiO
        moltype = 1
        we = 1241.5
        be = 0.7267
        ce = zero
        ae = 0.00494
        xmol = 1.0d-05
        popiso = 5.0d-02
!	NH3 v2 (average of s+a)
        moltype = 3
        we = 950.
        be = 9.947
        ce = 6.227
        ae = 0.0
!	for each of s, a
        xmol = 4.0d-06
!	a wild guess
        gcoll = 1.0d-13
        tcoll = 1.d+04
        popiso = 3.0d-03
!	HCN v2
        moltype = 2
        we = 712.
        be = 1.4872
        ce = zero
        ae = -8.0d-03
        xmol = 1.0d-06
        popiso = 1.6d-02
        acoeff = 2.43d-05
!	CO
        moltype = 1
        we = 2143.25
        be = 1.9225
        ce = zero
        ae = 0.01748
        xmol = 1.0d-04
        popiso = 1.6d-02
!	J=1-0 Einstein A
        acoeff = 7.2d-08

!	<sigma*v>/T^1/2
        grot = 1.0d-11
!	convert to 'band strength'
        arot = 3.*acoeff
        sv = arot/(8.*pi*clt*(two*be)**2)
!	emfac = two*hc2*w**3*sv*(xju/gju)*dju = (hc*w*acoeff/4*pi)*dju
!	ncrit = acoeff/gamma = (xju**4/gju)*arot/(grot*T**.5)

!	skip molecular lines if moltype = 0
!	moltype = 0

        wec = we/clt
        hcb = hc*be
        hcc = hc*ce
        hcwe = hc*we

!	default input parameters

!	outer radius of grid 1
        rgrid1 = 0.20*au
!	column density knee radius
        r0 = au
!	outer knee radius
        r1 = 400.*au
!	gas column density in disk at r0 (g/cm2 to be converted to mol/cm2)
        sigma0 = 100.0
!	column density of inner disk at r0
        sigma1 = 0.0
!	radial power law negative index
!	C&G use 1.5, but D''A finds ~1.0 for constant alpha
!	this is now for pres, so ~ density pcol + 1.5
        pcol = 2.5
!	radial power law positive index for column density inside r0
        qcol = 0.0
!	negative index for extra inner disk
        ocol = 2.5
!	to specify sigma at 1 AU instead of inside r0
!        sigma0 = sigma0*((au/r0)**pcol+(r0/au)**qcol)
!	turbulent viscosity parameter / dust density for dust settling
!	make negative to make temperature dependent
        avisc = 1.0d-02
        rhod = 3.0d+00
!	radius of star
        rstar = 2.0*rsun
!	mass of star (M_sun)
        xmstar = 1.0
!	luminosity of star (L_sun)
        xlstar = 4.0
!	add a hot source of UV to stellar spectrum
        tuv = 1.5d+04
        xuv = 0.1
!	carbon and silicon fraction in grains
        xc = 0.1
        xsi = one
!	carbon and silicate rolloff radius (cm) (amax = 0.1 cm)
        ac = 1.0d-05
        asi = 1.0d-05
!	dust sublimation temperature coefficients (for both carbon and silcate)
!	tsubli = tsubl*pres**gsubl
        tsubl = 600.
        gsubl = 0.02
!	destroy carbon grains if tgas > tburn
        tburn = 0.
!	enhance CO where carbon is burned
        xburn = 1.
!	enhancement of stellar radiation by reflection and emission from disk
        enhanc = 0.5
!	accretion rate for heating
        xmdot = 1.0d-08
!	dust to gas heat flow acommodation coefficient
        acom = 0.05d+00
!	x-ray luminosity of star
        pxr = 0.0d-04
!	x-ray absorption cross section per H
!	~1e-22 at 1 keV, falling as ~E^-2.5 (Gorti&Hollenbach)
        axr = 1.0d-22
!	extra x-ray parameter not currently used
        bxr = 2.5
!	photoelectric heating of gas
        npe = nu-4
        epspe = 0.05d+00
!	abundances for cooling function; should vary in disk
        xco0 = 1.0d-04
        xwat0 = 1.0d-04
        tfrco = 3.0d+01
        tfrwat = 1.5d+02
!	angular width of wind-swept cone
        coswind = 0.5
!	mass-loss rate in wind
        wmdot = 0.0d-10
!	wind speed
        vwind = 10.0
!	molecule destruction temperature
!	CO should extend in to where T = 4000.
        tdiss = 4000.
!	molecule freeze-out temperature
        tfrmol = 50.
!	destroy molecules where exposed to stellar UV or x-rays
!	this is now the radial column, not the vertical
!	GNI04 say 1e21 (vertical) for CO, and larger for bigger molecules
        coldiss = 1.0d+22
!	assume H is atomic, so increase collision rate, out to colHI
!	GNI04 say 4e21 (vertical), but PDRs have colCO ~ colHI
!	not currently used
        colhi = 1.0d+22
!	fluorescent excitation band strength
!	multiply by gguv = gg(iufl)
        sfluor = zero
        wfl = 5.0d+04
        fluor = fpi*sfluor/(hc*wfl)
!	microturbulent velocity parameters
        dvturb0 = 0.25*vkms
        dvturb1 = 0.25*vkms
!	velocity step
        dv = 0.05*vkms

!	assume LTE level populations?
        lte = .true.
!	C&G iterations
        kit = 4
!	dust radiative transfer iterations
        lit = 16
!	maximum line radiative transfer iterations (currently ignored)
        mit = 16
!	number of grain sizes or types; first half are carbon
!	parameter (ma = 32)
        na = 24
        nac = 12
!	number of wavelengths for dust spectrum
!	parameter (mu = 64)
        nu = 64
!	number of wavelengths for line spectra
!	parameter (mw = 1024)
        nw = mw
        nw2 = nw/2
        wmid = (nw+one)/two
!	number of nested grids (<= 12)
!	parameter (mg = 20)
        ng = 12
!	grid density multiplier for inner grid
!	nn9 = 8 make pretty pictures, but it takes forever for the lines
        nn9 = 4
!	parameter (mr = 128)
        nr = 128
!	costh runs from 1 to 0
!	set nnt = nr so r*dcosth = dr at outer edge of grids
!	or maybe bigger is better (and if mult of 3 can use 12 processors)
        nnt = 3*nr/2
!	parameter (mt = 192)
!	don''t bother calculating at theta above disk
        nt = nnt/2
!	need to calculate all nnt if wmdot /= 0
        if(wmdot.gt.zero) nt = nnt
!	nx = 16 seems to be OK for SED, but maybe not for disk structure
        nx = 64
!	use lx,lrot for nx,nz,nrot for the line calculation
        lx = 16
!	make nrot a fraction of nnt (or vice versa?) and a multiple of mysize
        nrot = nnt/2
        if(mod(nrot,mysize).ne.0) then
          if(doprint) print '(" nrot,mysize =",2i6)', nrot,mysize
          if(ifmpi) call mpi_abort(MPI_COMM_WORLD,icode,ierr)
          stop
        endif
        lrot = nrot/4
        smooth = .true.
        if(smooth) then
!	angular smoothing of radiation field in grid 9.  ntsm = 2 in grids 1-8
          ntsm = 2
!	radial smoothing hardwired for nrsm = 2 if set > 1
          nrsm = 2
        else
          ntsm = 1
          nrsm = 1
        endif
        smooth = .false.  !  turn off extra angular smoothing in grid 9
!	wavenumber (>0) or grain type (<0) for ggfs plots
        iupl = -6
!	outermost grid in cubes
        ngc = ng/2
!	grid for ggfs plots
        igpl = 6
!	parameter (mc = 256)
!	parameter (mmap = 1024)
        nc = mc

        call readparms

!	stellar effective temperature
        if(xlstar*rstar.gt.zero) then
          tstar = 5.8d+03*(xlstar/(rstar/rsun)**2)**quart
        elseif(xlstar*tstar.gt.zero) then
          rstar = rsun*sqrt(xlstar)*(5.8d+03/tstar)**2
        elseif(rstar*tstar.gt.zero) then
          xlstar = (tstar/5.8d+03)**4*(rstar/rsun)**2
        endif
        if(tburn.le.zero) then
          tburn = tstar
          xburn = one
        endif
!	T s.t. kT = G Mstar Mmol / Rstar
        tgrav = 6.67d-08*xmstar*xmsun*xmg/(rstar*bk)
!	effective temperature including UV excess
        teff = tstar*(one+xuv)**quart
!	hot dust temperature at rstar
        tsurf = teff/(two**0.4)
        gmmdot = 8.36d+51*xmstar*xmdot
        wmdot = wmdot*xmsun/(xmg*secyr*fpi)
        if(doprint.and.(wmdot.gt.zero).and.(nt.lt.nnt)) &
          print '("Warning: nt must = nnt if wmdot > 0")'
!	viscosity parameter for dust settling
        avisc = avisc*xmg/rhod
!	x-ray luminosity of star * heating efficiency
        pxr = 0.4*pxr*xlstar*xlsun
!	orbital speed at 1 AU (scale as r^-1/2)
        vorb = 30.*vkms*sqrt(xmstar)
!	vterm may be hardwired = 0
        vterm = zero*vkms
        iwterm = nw2+one+vterm/dv
!	wind speed
        vwind = vwind*vkms
!	Gaussian resolution HW1/e in pixels
        dsm = one*vkms/dv
        nsm = ten*dsm
        nsm = min(32,nsm)
        ngr = ng+nn9
        xn9 = nn9
        xnn9 = xn9*xn9
        ng9 = ng+1
        nr9 = nn9*nr
!	cell size on grid 0; dr = dr0*2**ig
        dr0 = rgrid1/(2*nr)
!	cell size on expanded grid 1
        dr9 = two*dr0/xn9
        nr2 = nr/2
        nt8 = nt-nt/nn9
        dcost = one/nnt
        dcos9 = dcost/xn9
        coswind = min(coswind,(nt*dcost))
        thwind = acos(coswind)
        dcrot = one/nrot
!	calculate lower half of grid and mirror
        bipolar = ((vterm.eq.zero).and.(wmdot.eq.zero))
        nx9 = nn9*nx
        ny = nx
!	set nz = 2*nr so dz = dr
        nz = 2*max(nr,nnt)
!	need greater nz with inner disk?
!        if(sigma1.gt.zero) nz = 2*nz
        nz9 = nn9*nz
        xnx = nx
        xnz = nz
!	rho-z grid size
        nrz = 4*nr
!	convert sigma0 (roughly) to pmid0
        sigma0 = sigma0/xmg
        sigma1 = sigma1/xmg
        tmid0 = 3.00d+02
        hh2 = (tmid0/tgrav)*(r0**3/rstar)
        pmid0 = tmid0*sigma0/sqrt(tpi*hh2)
        pmid1 = tmid0*sigma1/sqrt(tpi*hh2)
!	pmid = (pmid0/((r0/r)**qcol+(r/r0)**pcol)
!		+pmid1/(r0/r)**ocol)*sexp(-r/r1)
!	inner and outer edge of calculated disk
        pmin = one
        rin = rgrid1/1.6d+01
        if(qcol.gt.zero) rin = max(rin,r0*(pmin/pmid0)**(one/qcol))
        if(sigma1.gt.zero) rin = 1.001*rstar
        rout = rgrid1*two**(ng-1)
        if(r0*pcol.gt.zero) then
          rmax = rout
          do iit = 1,4
            rmax = rmax*slog((pmid0/pmin)*(r0/rmax)**pcol)
          enddo
          rout = min(rout,rmax)
        endif

        npl = lrot
        if(npl.gt.mu) npl = npl/2
        if((ngr.gt.mg).or.(nr.gt.mr).or.(nt.gt.mt).or.(nrot.gt.2*mt) &
          .or.(npl.gt.mu)) then
          if(doprint) then
            print '(" array out of bounds",10i6,2i8)', &
            ngr,mg,nr,mr,nt,mt,nrot,mt,nz,nt*ng*nr2,mmap
            if(ifmpi) call mpi_abort(MPI_COMM_WORLD,icode,ierr)
            stop
          endif
        endif
        if(ifmpi.and.doprint &
          .and.(mod(nrot,mysize).ne.0)) then
          print '(" nrot must be a multiple of mysize",2i8)', &
            nrot,mysize
          if(ifmpi) call mpi_abort(MPI_COMM_WORLD,icode,ierr)
          stop
        endif

!	pass parameters to startgr
        call initgr(dr0,dr9,dcost,dcos9,nr,nr9,ng9,nt,nt8)

        dx0 = (two*nr*dr0)/xnx
        dx9 = two*dx0/xn9
        dz0 = xnx*dx0/xnz
        dz9 = two*dz0/xn9
        nx2 = nx/2
        ny2 = ny/2
        nz2 = nz/2
        if(doprint) then
!	typical values
          tgi = 300.
          hh2 = (tgi/tgrav)*(r0**3/rstar)
          dens = sigma0/sqrt(tpi*hh2)
          print '(" H2 dens, mass column, scale height at",f8.3, &
            " AU =",1p,3e9.2)', &
            r0/au,dens,sigma0*xmg,sqrt(hh2)/au
          print '(" stellar M, L, T, R, dMdt =", &
            2f8.3,f8.0,f8.2,1p,e9.2)', &
            xmstar,xlstar,tstar,rstar/rsun,xmdot
          print '(" vorb, vturb, vstep =",f8.1,3f8.3)', &
            vorb/vkms,dvturb0/vkms,dvturb1/vkms,dv/vkms
          print '(" inner xyz grid:",i3," x",f10.5," to",f10.4," AU")', &
            nx9,dx9/au,nx9*dx9/(two*au)
          print '(" outer xyz grid:",i3," x",f10.4," to",f10.3," AU")', &
            nx,(two**ng)*dx0/au,(two**(ng-1))*nx*dx0/au
          print '(" calculation limits (AU):",2f10.3)', &
            rin/au,rout/au
          write(16,'(" H2 dens, mass column, scale height at",f8.3, &
            " AU =",1p,3e9.2)') &
            r0/au,dens,sigma0*xmg,sqrt(hh2)/au
          write(16,'(" stellar M, L, T, R, dMdt =", &
            2f8.3,f8.0,f8.2,1p,e9.2)') &
            xmstar,xlstar,tstar,rstar/rsun,xmdot
          write(16,'(" vorb, vturb, vstep =",f8.1,3f8.3)') &
            vorb/vkms,dvturb0/vkms,dvturb1/vkms,dv/vkms
          write(16,'(" inner xyz grid:",i3," x",f10.5, &
            " to",f10.4," AU")') &
            nx9,dx9/au,nx9*dx9/(two*au)
          write(16,'(" outer xyz grid:",i3," x",f10.4, &
            " to",f10.3," AU")') &
            nx,(two**ng)*dx0/au,(two**(ng-1))*nx*dx0/au
        endif
        astar = pi*(rstar/(dx0*two**ng))**2
        dxstar = rstar/four
        gsum = zero
        do ix = 1,4
          x = dxstar*(ix-half)
          gwt = sqrt(one-(x/rstar)**2)
          gsum = gsum+two*gwt
        enddo
        bstar = astar/gsum

!	set up dust arrays
        tsubli = tsubl*max(one,pmid0)**gsubl
        call dustinit(xc,xsi,ac,asi,tuv,xuv,emiss,tsubli)
!	abss is the fraction of the incident starlight absorbed by the disk
        abss = sqrt(emiss)
!	replace bus with hd100 spectrum, normalized to same power
        sumbus = zero
        sumdat = zero
        do iu = 1,nu
          sumbus = sumbus+bus(iu)*ud(iu)
          sumdat = sumdat+spec(iu,15)
        enddo
        sumrat = sumbus/sumdat
        fwmx = zero
        sumfs = zero
        do iu = 1,nu
          bus(iu) = sumrat*spec(iu,15)/ud(iu)
          busdud = bus(iu)*dud(iu)
          if(busdud.gt.fwmx) fwmx = busdud
          sumfs = sumfs+busdud
          gg(iu,1,1,1) = zero
          opa(iu,1,1,1) = zero
        enddo
        fwmx = astar*fwmx
        if(doprint) then
          xlvis = sigsb*tstar**4*fpi*rstar**2/xlsun
          xluv = xuv*xlvis
          xlphot = sumfs*pi*fpi*rstar**2/xlsun
          xlaccr = gmmdot/(two*rin*xlsun)
          print '(" Lvis,Luv,Lphot,Laccr(rin,rstar) =",5f8.3)', &
            xlvis,xluv,xlphot,xlaccr,xlaccr*(rin/rstar)
!	heating and cooling power per atom at r0
          colin = 1.0d+22
          colup = 1.0d+21
          xco = xco0
          xwat = xwat0
          cgas = cooling(dens,tgi,1,1,1)
          popsum = zero
          do ia = 1,na
            popa = popd(ia)*pi*ad(ia)**2
            popsum = popsum+popa
          enddo
          acgd = four*dens**2*popsum*sqrt(bk*tgi/xmg)*acom*2.5*bk
          hxr = (pxr/(fpi*r0**2))*((axr/two)/(axr*colin/two+one)**2)
          haccr = gmmdot/(fpi*r0**3*sigma0)
          print '(" gas heating and cooling at r0:",1p,4e10.2)', &
            hxr,haccr,cgas,acgd
        endif
        if(reuse) then
          lw = nu
          if(doprint) print '(" reading saved disk model")'
          call readmodel
          if(doplot) then
            print '(" plotting log_10(gas density)")'
            call plotcon(dg,.true.,'')
            call waitasec(wait,2)
            print '(" plotting log_10(gas temperature)")'
            call plotcon(tg,.true.,'')
            call waitasec(wait,2)
          endif
          goto 300
        endif

        if(ifmpi) call mpi_barrier(MPI_COMM_WORLD,ierr)
        if(doprint) print '(" calculate C&G disk model")'

!	initialize taufac
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            do it = 1,nt
              opl(ml,it,ir,ig) = one
            enddo
          enddo
        enddo
!	calculate sublimation front and C&G disk structure
!	don''t know how to farm out
        ll = 1
        lw = 1
        ffstar = one
        nit = 0
!	rim location on the midplane
!	approximations may break down if qcol is too big
        powz = gsubl/0.2  !  0.1
        if(pmid1.gt.zero) then
          ginv = 0.4/gsubl+ocol  !  20+ocol
          powr = two+ocol*powz  !  2+ocol/10
        elseif(qcol.lt.ten) then
          ginv = 0.4/gsubl-qcol  !  20-qcol
          powr = two-qcol*powz  !  2-qcol/10
        else
          ginv = 0.4/gsubl-ten  !  10
          powr = two-ten*powz  !  1
        endif
        ustar = four*tstar/hck
        iustar = slog(ustar/ud(1))/slog(ud(2)/ud(1))+1.5
        iustar = min(nu,max(1,iustar))
        ias1 = na/2+1
        ias = slog(asi/ad(ias1))/slog(ad(ias1+1)/ad(ias1))+1.5
        ias = min(na,max(ias1,ias))
        r = rin
        pmid = pmid0/((r0/r)**qcol+(r/r0)**pcol)+pmid1/(r/r0)**ocol
        if(r1.gt.zero) pmid = pmid/sexp(r/r1)
        tsubli = tsubl*max(one,pmid)**gsubl
        usubl = four*tsubli/hck
        iusubl = slog(usubl/ud(1))/slog(ud(2)/ud(1))+1.5
        iusubl = min(nu,max(1,iusubl))
        rim0 = (rstar/two)*(teff/tsubli)**2 &
          *(phab(iustar,ias)/phab(iusubl,ias))**(one/powr)
        if(rim0.lt.rin) then
          if(doprint) print '("Warning: rim0 inside of rin")'
          rim0 = rin
        endif
        rim = rim0*(one+enhanc)**(one/powr)
!	calculate sublimation rim in wind
        if(wmdot*vwind.gt.zero) then
          pwind = tsubli*wmdot/(vwind*rim0**2)
          tsubli = tsubl*max(one,pwind)**gsubl
          usubl = four*tsubli/hck
          iusubl = slog(usubl/ud(1))/slog(ud(2)/ud(1))+1.5
          iusubl = min(nu,max(1,iusubl))
          poww = two+two*powz
          rimw = (rstar/two)*(teff/tsubli)**2 &
            *(phab(iustar,ias)/phab(iusubl,ias))**(one/poww)
!	calculate midplane wind-disk pressure balance rim
          rimb = ((wmdot*xmg*vwind)/(pmid0*bk))**(one/(qcol+two)) &
            *r0**(qcol/(qcol+two))
        else
          rimw = rout
          rimb = zero
        endif
        if(doprint) &
          print '(" initial disk, wind rims:",3f10.2)', &
          rim/au,rimw/au,rimb/au

        grmax = 0.99
!	don''t need lastit in frisco since don''t recalculate SED
        xmdisk = zero
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          if(igr.le.1) then
            ifrim = .false.
            zr1 = zero
            gr1 = one
            graz = one
            zrim = zero
            rphot = zero
            zphot = zero
            zrlast = zero
            do it = 1,nt
              do iu = 1,nu
                taud(iu,it) = zero
              enddo
              dtac4(it) = zero
            enddo
            itwind = nt
            itrim0 = nt
            itrim1 = nt
          endif
          do ir = kr,lr
            r = dr*(ir-half)
            if((r.gt.rin).and.(r.lt.rout)) then
              tst4 = ostar*teff**4
              tir4 = graz*abss*tst4/two
              tir = tir4**quart
!	use tir (leaving out dtac4) for calculation of sigma and graz?
!	and just use first term in pmid for sigma so heating doesn't blow up
              hr2 = (tir/tgrav)*(r/rstar)
              hh2 = hr2*r**2
              sigma = (pmid0*(r0/r)**pcol/tir)*sqrt(tpi*hh2)
              sigma = max(sigma,1.0d+18)
              pmid = pmid0/((r0/r)**qcol+(r/r0)**pcol) &
                     +pmid1/(r/r0)**ocol
              if(r1.gt.zero) pmid = pmid/sexp(r/r1)
              tsubli = tsubl*max(one,pmid)**gsubl
            else
              sigma = zero
              pmid = zero
              tsubli = tsubl
            endif
!	find z rim semianalytically
            if(r.lt.rim) then
              if((pmid.gt.one).and.(ig.ne.1)) then
                usubl = four*tsubli/hck
                iusubl = slog(usubl/ud(1))/slog(ud(2)/ud(1))+1.5
                iusubl = min(nu,max(1,iusubl))
                rim0 = (rstar/two)*(teff/tsubli)**2 &
                  *(phab(iustar,ias)/phab(iusubl,ias))**(one/powr)
                rim = rim0*(one+enhanc)**(one/powr)
              endif
              zri = zero
              gri = one
              prim = pmid
            else
              if(doprint.and.(.not.ifrim)) print &
                '(" found sublimation rim at",3f8.4,f8.2,1p,e10.2)', &
                rim/au,rim0/au,rin/au,tsubli,pmid/tsubli
              ifrim = .true.
              alogr = slog(r/rim0)
              rr0 = (r/rim0)**powr
!	do I want to use gr1, or graz, here?
              arg = max(zero,ginv*alogr-slog(one+gr1*enhanc)/powz)
              if(zr1.eq.zero) then
!	first point beyond rim
                zri = sqrt(arg*hr2)
              else
!	z rim with last value of gr1
                zri = sqrt(arg*hr2)
!	solve for z rim iteratively
                do iit = 1,4
!	sin(graz)
                  gri = (rr0*sexp(-powz*zri**2/hr2)-one)/enhanc
                  gri = max(zero,min(grmax,gri))
!	tan(graz)
                  gri = gri/sqrt(one-gri**2)
                  zri = (zri+zr1+gri*dr/r)/two
                enddo
!	z rim including disk radiation
                gri = gri/sqrt(one+gri**2)
                arg = ginv*alogr-slog(one+gri*enhanc)/powz
                if(arg.gt.zero) zri = sqrt(arg*hr2)
              endif
              zrim = zri*r
              if(arg.gt.zero) then
                prim = pmid*exp(-arg)
              else
                prim = pmid
              endif
            endif
            zr1 = zri
            gr1 = gri
            zphot = zero
            itrim0 = itrim1
            do it = 1,nt
              costh = dcosth*(nt-it+half)
              z = r*costh
              if((rphot.gt.zero).and.(r.gt.rphot).and.(zphot.eq.zero) &
                .and.(taud(iuv,it).gt.one)) then
                if((ig.eq.ng9).and.(it.eq.1)) then
                  zphot = r*min(coswind,zrlast*sqrt(slog(r/rphot)))
                elseif((ig.lt.ng9).and.(costh.ge.coswind)) then
                  if(doprint) print '("Warning: dust opaque at thwind", &
                    3i6)', it,ir,ig
                  zphot = r*coswind
                else
                  dit = (taud(iuv,it)-one)/(taud(iuv,it)-taud(iuv,it-1))
                  zphot = r*dcosth*(nt-it+half+dit)
                  zrlast = (zphot/r)/sqrt(slog(r/rphot))
                endif
              endif
!	initialize condens = 1
!	will reset condens = 0 at p < prim
!	burn carbon if Tg > tburn only after calculating temperature
              do ia = 1,na
                condens(ia,it,ir,ig) = one
              enddo
            enddo  !  it
            if((rphot.eq.zero).or.(r.lt.rphot).or.(zphot.eq.zero)) then
              graz = one
            else
              graz = zphot/(two*r*slog(r/rphot))
              graz = graz/sqrt(one+graz**2)
            endif
            if(r.le.rstar) then
              r = rstar
              ostar = one
            elseif(r.lt.(hun*rstar)) then
              ostar = (one-sqrt(one-(rstar/r)**2))/two
            else
              ostar = quart*(rstar/r)**2
            endif
            tdrad = teff*ostar**fifth
            tst4 = ostar*teff**4
            tir4 = graz*abss*tst4/two
            tir = tir4**quart
            hr2 = (tir/tgrav)*(r/rstar)
            hh2 = hr2*r**2

!	C&G interior temperature and density
!	wind density at r
            if(wmdot*vwind.ne.zero) then
              dgw = wmdot/(vwind*r**2)
            else
              dgw = zero
            endif
!	wind ram pressure
            pram = dgw*xmg*vwind**2/bk
!	iterate at each r
            nit = 0
  100       nit = nit-1
            pres = pmid
            ddfac = zero
            faccr = zero
            inwind = (pmid.le.pram).and.(itwind.eq.nt)
!	step up from midplane
            do it = nt,1,-1
              costh = dcosth*(nt-it+half)
              sinth = sqrt(one-costh**2)
              if(costh.gt.coswind) inwind = .true.
              z = r*costh
              dz = r*dcosth
  110         if((ig.eq.1).and.(it.gt.nt8)) then
!	average arrays from grid 9
                jr2 = ir*nn9
                jr1 = jr2-nn9+1
                jt2 = nt-(nt-it)*nn9
                jt1 = jt2-nn9+1
                sumdg = zero
                sumdd = zero
                sumdm = zero
                sumtg = zero
                sumtd = zero
                sumph = zero
                sumopl = zero
                sumcol = zero
                do jt = jt1,jt2
                  do jr = jr1,jr2
                    sumdg = sumdg+dg(jt,jr,ng9)
                    sumdd = sumdd+dd(jt,jr,ng9)
                    sumdm = sumdm+dm(jt,jr,ng9)
                    sumtg = sumtg+tg(jt,jr,ng9)
                    sumtd = sumtd+td(jt,jr,ng9)
                  enddo
                  sumopl = sumopl+opl(ml,jt,jr2,ng9)
                  sumcol = sumcol+colr(jt,jr2,ng9)
                enddo
                dens = sumdg/xnn9
                ddfac = sumdd/xnn9
                tgi = sumtg/xnn9
                tdi = sumtd/xnn9
                pres = dens*tgi
                iamax = int(sumdm/xnn9)
                opl(ml,it,ir,ig) = sumopl/xn9
                colr(it,ir,ig) = sumcol/xn9
                if(sigma.gt.zero) then
                  faccr = faccr+gmmdot*dens*dz/(fpi*sigma*r**3)
                endif
                opapl = zero
                bwsum = zero
                do ia = 1,na
                  sumcon = zero
                  do jt = jt1,jt2
                    do jr = jr1,jr2
                      sumcon = sumcon+condens(ia,jt,jr,ng9)
                    enddo
                  enddo
                  condens(ia,it,ir,ig) = sumcon/xnn9
                enddo
                do iu = 1,nu
                  sumopa = zero
                  sumsca = zero
                  do jt = jt1,jt2
                    do jr = jr1,jr2
                      sumopa = sumopa+opa(iu,jt,jr,ng9)
                      sumsca = sumsca+sca(iu,jt,jr,ng9)
                    enddo
                  enddo
                  opai = sumopa/xnn9
                  opa(iu,it,ir,ig) = opai
                  scai = sumsca/xnn9
                  sca(iu,it,ir,ig) = scai
                  if(opai.gt.zero) then
                    bwdud = bw(ud(iu),tgi)*dud(iu)
!	since bw/opa ~ w, not w**3
                    if(iu.eq.1) bwdud = two*bwdud
                    opapl = opapl+bwdud/opai
                    bwsum = bwsum+bwdud
                  else
                    opapl = zero
                    goto 112
                  endif
                enddo
                opapl = bwsum/opapl
                if(badnum(opapl)) then
                  if(doprint) print '("at",3i6," opapl =",1p,5e10.2)', &
                    it,ir,ig,opapl,bwsum,tgi,bwdud,opai
                  opapl = zero
                endif
  112           dtac4(it) = thrq*faccr*opapl*dz/sigsb
              else  !  (ig.ne.1).or.(it.le.nt8)
!	start with temperature from last iteration
!	or extrapolated from neighboring cells
                if(nit.lt.(-1)) then
                  tgi = tg(it,ir,ig)
                  tdi = td(it,ir,ig)
                elseif(ir.gt.kr) then
                  if((it.eq.nt).or.((ig.eq.1).and.(it.eq.nt8))) then
                    tgi = tg(it,ir-1,ig)
                    tdi = td(it,ir-1,ig)
                  elseif((it.le.itrim0).and.(it.ge.itrim1)) then
                    tgi = tg(it+1,ir,ig)
                    tdi = td(it+1,ir,ig)
                  else
                    tgj = tg(it+1,ir,ig)
                    tgk = tg(it,ir-1,ig)
                    tgl = tg(it+1,ir-1,ig)
                    tgi = max(tgj,(tgj+tgk-tgl))
                    tdj = td(it+1,ir,ig)
                    tdk = td(it,ir-1,ig)
                    tdl = td(it+1,ir-1,ig)
                    tdi = max(tdj,(tdj+tdk-tdl))
                  endif
                elseif(igr.gt.1) then
                  if(it.eq.nt) then
                    tgi = tg(it,lr,ig-1)
                    tdi = td(it,lr,ig-1)
                  elseif((it.le.itrim0).and.(it.ge.itrim1)) then
                    tgi = tg(it+1,ir,ig)
                    tdi = td(it+1,ir,ig)
                  else
                    tgj = tg(it+1,ir,ig)
                    tgk = (tg(it,lr-1,ig-1)+tg(it,lr,ig-1))/two
                    tgl = (tg(it+1,lr-1,ig-1)+tg(it+1,lr,ig-1))/two
                    tgi = max(tgj,(tgj+tgk-tgl))
                    tdj = td(it+1,ir,ig)
                    tdk = (td(it,lr-1,ig-1)+td(it,lr,ig-1))/two
                    tdl = (td(it+1,lr-1,ig-1)+td(it+1,lr,ig-1))/two
                    tdi = max(tdj,(tdj+tdk-tdl))
                  endif
                else
                  tgi = tstar
                  tdi = tstar
                endif
                tgi0 = tgi
                if(r.lt.rin) then
                  dens = zero
                  pres = zero
                  ddfac = zero
                  tsubli = tsubl
                  tdhot = tsubli/two**quart
                  do iu = 1,nu
                    trans = one
                    gg(iu,it,ir,ig) = &
                      bw(ud(iu),tdhot)*(tir4/tdhot**4)*trans/two &
                      +ostar*(bus(iu)+graz*rebus(iu))*trans
                  enddo
                elseif(inwind) then
                  dens = dgw
                  pres = dens*tgi
                  ddfac = zero
                  tsubli = tsubl*max(pres,one)**gsubl
                  tdhot = tsubli/two**quart
                  do iu = 1,nu
                    trans = sexp(-taud(iu,it))
                    gg(iu,it,ir,ig) = &
                      bw(ud(iu),tdhot)*(tir4/tdhot**4)*trans/two &
                      +ostar*(bus(iu)+graz*rebus(iu))*trans
                  enddo
                else
                  hr2 = (tgi/tgrav)*(r/rstar)
                  hh2 = hr2*r**2
!	this may work.  I should try it.
!                  if(badnum(pgj).or.badnum(pgk).or.badnum(pgl)) &
!                    print '(3i6,1p,3e10.2)', it,ir,ig,pgj,pgk,pgl
!                  pgj = slog(max(one,pgj))
!                  pgk = slog(max(one,pgk))
!                  pgl = slog(max(one,pgl))
!                  dth = acos(costh-dcosth)-acos(costh)
!                  rdth = (r-dr)*dth
!                  dpt = (pgk-pgl)/rdth
!                  dpr = (pgj-pgl)/dr
!                  dprh = sinth*dpr-costh*dpt
!                  dpz = -z/hh2
!                  dpr = (costh*dpz+sinth*dprh)*dr
!                  dpt = (sinth*dpz-costh*dprh)*rdth
!	should make this change more gradually
!                  if(costh.gt.half) then
!                    pres = sexp(pgk+dpr)
!                    dpt = pgk+dpr-pgj
!                  else
!                    pres = sexp(pgj+dpt)
!                  endif
!	pressure change in half cell
                  gaus = sexp(-(z*dz)/(two*hh2))
!	pressure at center of cell
                  pres = pres*gaus
                  if(badnum(pres).or.(pres.lt.one)) pres = zero
                  condi = condens(na/4,it,ir,ig)
                  if(pres.le.prim) then
                    do ia = 1,na
                      condens(ia,it,ir,ig) = zero
                    enddo
                  else
                    if(it.le.itrim1) itrim1 = it-1
                    do ia = 1,na
                      condens(ia,it,ir,ig) = one
                    enddo
                  endif
                  if((it.le.itwind).and.(pres.lt.pram)) then
                    inwind = .true.
                    itwind = it
                    goto 110
                  endif
!	density at center of cell
                  dens = pres/tgi
!	dust density decrease exponential factor for a = 1 cm
!	(avisc has been multiplied by xmg/rhod)
                  if((avisc.gt.zero).and.(dens.gt.zero)) then
                    ddfac = ddfac+z*dz/(avisc*hh2*sqrt(hh2)*dens)
                  elseif((avisc.lt.zero).and.(dens.gt.zero)) then
                    avist = -avisc*min(ten,max(one,tgi/hun))**4
                    ddfac = ddfac+z*dz/(avist*hh2*sqrt(hh2)*dens)
                  else
                    ddfac = zero
                  endif
!	accumulated heat flux from viscous heating
                  if(sigma.gt.zero) then
                    faccr = faccr+gmmdot*dens*dz/(fpi*sigma*r**3)
                  endif
!	add accretion heating to heating by hot dust emission
!	tir4 is from midplane calculation
!	dtac4 is from previous r or previous iteration
                  tir = (tir4+dtac4(it))**quart
                  tsubli = tsubl*max(pres,one)**gsubl
                  tdhot = tsubli/two**quart
                  tir = min(tir,tsubli)
                  do iu = 1,nu
                    transr = sexp(-taud(iu,it))
                    trans = sexp(-taud(iu,it)*graz)
!	does this include downward scattered light?
                    gg(iu,it,ir,ig) = bw(ud(iu),tir)*(one-trans/two) &
                      +bw(ud(iu),tdhot)*(tir4/tdhot**4)*trans/two &
                      +ostar*(bus(iu)*transr+graz*rebus(iu)*trans)
                  enddo
                endif  !  inwind
                taucont = taud(iuv,it)
                dg(it,ir,ig) = dens
                dd(it,ir,ig) = ddfac
                if(ir.gt.kr) then
                  colin = colr(it,ir-1,ig)
                  colup = colz(it,ir-1,ig)
                elseif(igr.gt.1) then
                  colin = colr(it,lr,ig-1)
                  colup = colz(it,lr,ig-1)
                else
                  colin = zero
                  colup = zero
                endif
!	should use better chemistry
                xwat = xwat0*sexp(-(tfrwat/tdi)**10)
                xco = xco0*sexp(-(tfrco/tdi)**10)
!	enhance CO where carbon is burned
                if((xburn.gt.one).and.(tburn.gt.one)) &
                  xco = xco*(xburn+(one-xburn)*sexp(-(tgi/tburn)**10))
                if(taucont.le.zero) then
                  xwat = xwat/hun
                  xco = xco/hun
                endif
!	calculate dust temperatures and opacities
                call gastemp(it,ir,ig,iamax,tsubli,tgi,tdi,hcgas)
                dens = pres/tgi
                if(dens.lt.tiny) then
                  dens = zero
                  pres = zero
                endif
!	calculate Planck mean opacity (should use Rosseland?)
                opapl = zero
                bwsum = zero
                do iu = 1,nu
                  opai = opa(iu,it,ir,ig)
                  if(opai.gt.zero) then
                    bwdud = bw(ud(iu),tgi)*dud(iu)
                    if(iu.eq.1) bwdud = two*bwdud
                    opapl = opapl+bwdud/opai
                    bwsum = bwsum+bwdud
                  else
                    opapl = zero
                    goto 120
                  endif
                enddo
                opapl = bwsum/opapl
                if(badnum(opapl)) then
                  if(doprint) print '("at",3i6," opapl =",1p,5e10.2)', &
                    it,ir,ig,opapl,bwsum,tgi,bwdud,opai
                  opapl = zero
                endif
  120           dtac4(it) = thrq*faccr*opapl*dz/sigsb
                taufac = opl(ml,it,ir,ig)
                if(.not.inwind) then
                  hr2 = (tgi/tgrav)*(r/rstar)
                  hh2 = hr2*r**2
!	undo first guess at pressure change
                  pres = pres/gaus
!	pressure change in full cell
                  gaus = sexp(-(z*dz)/hh2)
!	pressure at top of cell
                  pres = pres*gaus
                  if(pres.lt.one) pres = zero
                endif
              endif  !  ig.eq.1?
!	values at cell center
              dg(it,ir,ig) = dens
              dd(it,ir,ig) = ddfac
              tg(it,ir,ig) = tgi
              td(it,ir,ig) = tdi
              colr(it,ir,ig) = colin+taufac*dens*dr
              sumcon = zero
              do ia = 1,na
                sumcon = sumcon+condens(ia,it,ir,ig) &
                  *sexp(-ad(ia)*ddfac)*taufac
              enddo
              dm(it,ir,ig) = sumcon
!              dm(it,ir,ig) = iamax
!              dm(it,ir,ig) = tdi/tsubli
!              dm(it,ir,ig) = tdi/tgi
              if(pres.lt.one) dm(it,ir,ig) = one
!              dm(it,ir,ig) = sqrt(opa(iusil,it,ir,ig) &
!                *bw(ud(iusil),tdi)*r**3)
            enddo  !  it
!	save tmid for next iteration
            tmid = (9.*tg(nt,ir,ig)-tg(nt-1,ir,ig))/8.
!	integrate dtac4 and colz down through the disk
            dtaci = faccr/(two*sigsb)
            dtac4(1) = dtaci
            colup = zero
!	hh2 and z should still be known from last it (1)
            if(ig.eq.ng9) colup = dg(1,ir,ig)*(hh2/z)
            colz(1,ir,ig) = colup
            do it = 2,nt
              dtaci = dtaci+dtac4(it)
              dtac4(it) = dtaci
!	leave out opl(ml) since stretching is radial
              colup = colup+dg(it,ir,ig)*dz
              colz(it,ir,ig) = colup
            enddo
            if(nit.gt.-kit) goto 100

            xmdisk = xmdisk+fpi*r*dr*colup
            iw = ir+(ig-1)*nr/2
            if((ig.le.ng).and.(iw.ge.1).and.(iw.le.mw)) then
              smax = 1.0d+04
              spec(iw,4) = slog10(max(one,min(smax,tg(nt,ir,ig))))
              spec(iw,16) = spec(iw,4)
            endif
            do it = 1,nt
              do iu = 1,nu
                opai = opa(iu,it,ir,ig)
                if(opai.gt.zero) then
                  taud(iu,it) = taud(iu,it) &
                    +sqrt(opai*(opai+sca(iu,it,ir,ig)))*dr
                endif
              enddo
            enddo
            if((rphot.eq.zero).and.(taud(iuv,nt).gt.one)) then
              dtau = (opa(iuv,nt,ir,ig)+sca(iuv,nt,ir,ig))*dr
              dir = (taud(iuv,nt)-one)/dtau
              rphot = dr*(ir+half-dir)
            endif
          enddo  !  ir
          if(ig.lt.ng) then
            tg(nt,nr/2,ig+1) = (tg(nt,nr-1,ig)+tg(nt,nr,ig))/two
            td(nt,nr/2,ig+1) = (td(nt,nr-1,ig)+td(nt,nr,ig))/two
          endif
        enddo  !  igr
        call copyarr(dm)
        call copyarr(dg)
        call copyarr(dd)
        call copywta(tg,dg)
        call copywta(td,dg)
        call copyarrr(opl,ml,ml)
        call copyarrr(opa,mu,nu)
        call copyarrr(sca,mu,nu)
        call copyarrr(swd,mu,nu)
        call copyarrr(gg,mu,nu)

        if(doprint) then
          print '(" disk mass:",1p,e10.2)', xmdisk*xmg/xmsun
          print '(" midplane sublimation, photosphere radii:", &
            2f10.4)', rim/au,rphot/au
          if(rphot.gt.rgrid1) &
            print '("Warning: photosphere outside of grid 1")'
          write(16,'(" midplane sublimation, photosphere radii:", &
            2f10.4)') rim/au,rphot/au
        endif

!	plot mean intensity / stellar on xy grid
        if(doplot) then
          print '(" plotting settling*condensation*taufac")'
!          print '(" plotting largest sublimating grain")'
!          print '(" plotting dust temperature / sublimation")'
!          print '(" plotting dust/gas temperature")'
!          print '(" plotting silicate emission * r^3")'
          call plotcon(dm,.false.,'')
          call waitasec(wait,2)
!	plot dust heating image
!          print '(" plotting dust heating radiation / stellar")'
!          if(rim.lt.rgrid1) then
!            call plotggfs(gg,ng9,iupl)
!            call waitasec(wait,2)
!          endif
!          call plotggfs(gg,igpl,iupl)
!          call waitasec(wait,2)
!	plot gas density and temperature
          print '(" plotting log_10(hydrogen density)")'
          call plotcon(dg,.true.,'')
          call waitasec(wait,2)
          print '(" plotting C&G gas temperature")'
          call plotcon(tg,.false.,'')
        endif

        if(ifmpi) call mpi_barrier(MPI_COMM_WORLD,ierr)
        lw = nu
        sumfs = sumfs*astar
!          +gmmdot/(two*max(rim,r0)*fpi*(dx0*two**ng)**2)

        nit = 0
        moltype = -abs(moltype)
        ffstar = one
        hcsum = zero

!	loopback point for dust radiative transfer
  200   nit = nit+1

        if(doprint) print '(" calculate hydro equil disk structure")'
!	calculate hydro equil disk structure with new temperatures
!	don''t know how to farm out
        rphot = zero
        xmdisk = zero
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          drho = dr/four
          dz = dr/four
          lp = 4*lr
          kp = lp/2+1
          lz = 4*nr
          kz = lz/2+1
          drig = dr
          lrig = lr
          rmx = (lr-half)*dr
          if(igr.le.1) cosw = zero
          do ip = 1,lp
            rho = (ip-half)*drho
!	calculate midplane temp, dens, pres
            r = rho
            tmid = (9.*tg(nt,ip,ig)-tg(nt-1,ip,ig))/8.
            hh2 = (tmid/tgrav)*(r**3/rstar)
            if((r.gt.rin).and.(r.lt.rout)) then
              pmid = pmid0/((r0/r)**qcol+(r/r0)**pcol) &
                     +pmid1/(r/r0)**ocol
              if(r1.gt.zero) pmid = pmid/sexp(r/r1)
            else
              pmid = zero
            endif
            dgmid = pmid/tmid
            if(badnum(dgmid)) dgmid = zero
            if(wmdot*vwind.ne.zero) then
              dgw = wmdot/(vwind*r**2)
            else
              dgw = zero
            endif
            pram = dgw*xmg*vwind**2/bk
            inwind = (pram.gt.pmid).and.(cosw.eq.zero)
!	solve eqn of hydro equil in cylindrical coords
            pres = pmid
            ddfac = zero
            colup = zero
!	use nr here b/c nz9 = nr, not nr9
            do iz = 1,lz
              z = (iz-half)*dz
              if((ig.eq.1).and.(iz.le.(lz/nn9))) then
!	interpolate from grid 9 onto grid 1
                jg = ng9
                jp2 = ip*nn9
                jp1 = jp2-nn9+1
                jz2 = iz*nn9
                jz1 = jz2-nn9+1
                sumnpz = xnn9
                sumtg = zero
                sumdg = zero
                sumdd = zero
                do jp = jp1,jp2
                  do jz = jz1,jz2
                    sumtg = sumtg+tgz(jz,jp,jg)
                    sumdg = sumdg+dgz(jz,jp,jg)
                    sumdd = sumdd+ddz(jz,jp,jg)
                  enddo
                enddo
                tgi = sumtg/sumnpz
                dens = sumdg/sumnpz
                ddfac = sumdd/sumnpz
!	want pres at top of cell
                jz1 = jz2-1
                jp1 = (jp1+jp2)/2
                jp2 = jp1+1
                pres1 = (dgz(jz1,jp1,jg)+dgz(jz1,jp2,jg)) &
                       *(tgz(jz1,jp1,jg)+tgz(jz1,jp2,jg))/four
                pres2 = (dgz(jz2,jp1,jg)+dgz(jz2,jp2,jg)) &
                       *(tgz(jz2,jp1,jg)+tgz(jz2,jp2,jg))/four
                pres = (3.*pres2-pres1)/two
              elseif((igr.gt.1).and.(ip.lt.kp).and.(iz.lt.kz)) then
!	interpolate from grid ig-1
                jg = ig-1
                jp2 = 2*ip
                jp1 = jp2-1
                jz2 = 2*iz
                jz1 = jz2-1
                sumnpz = four
                sumtg = zero
                sumdg = zero
                sumdd = zero
                do jp = jp1,jp2
                  do jz = jz1,jz2
                    sumtg = sumtg+tgz(jz,jp,jg)
                    sumdg = sumdg+dgz(jz,jp,jg)
                    sumdd = sumdd+ddz(jz,jp,jg)
                  enddo
                enddo
                tgi = sumtg/sumnpz
                dens = sumdg/sumnpz
                ddfac = sumdd/sumnpz
                pres1 = (dgz(jz1,jp1,jg)+dgz(jz1,jp2,jg)) &
                       *(tgz(jz1,jp1,jg)+tgz(jz1,jp2,jg))/four
                pres2 = (dgz(jz2,jp1,jg)+dgz(jz2,jp2,jg)) &
                       *(tgz(jz2,jp1,jg)+tgz(jz2,jp2,jg))/four
                pres = (3.*pres2-pres1)/two
              else
                r = sqrt(rho**2+z**2)
                costh = z/r
                xit = nt-costh/dcosth+half
                if(r.gt.rmx) then
!	can''t be more than one grid out from ig
                  if(igr.eq.0) then
                    jg = 2
                    dr = four*dr0
                    lr = nr
                    xit = nt-costh/dcost+half
                  elseif(igr.lt.ng) then
                    jg = igr+1
                    lr = lrig
                    dr = two*drig
                  else
                    jg = ig
                    lr = lrig
                    dr = drig
                  endif
                elseif(xit.lt.half) then
!	assume we''re in grid 9
                  jg = 1
                  dr = two*dr0
                  lr = nr
                  xit = nt-costh/dcost+half
                else
                  jg = ig
                  lr = lrig
                  dr = drig
                endif
                xir = r/dr+half
                ir = xir
                ir = max(1,min(lr-1,ir))
                jr = ir+1
                dir = xir-ir
!	don''t extrapolate
                dir = max(zero,min(one,dir))
                djr = one-dir
                it = xit
                it = max(1,min(nt-1,it))
                jt = it+1
                dit = xit-it
                dit = max(zero,min(one,dit))
                djt = one-dit
                tgi = djr*(djt*tg(it,ir,jg)+dit*tg(jt,ir,jg)) &
                      +dir*(djt*tg(it,jr,jg)+dit*tg(jt,jr,jg))
                taufac = djr*(djt*opl(ml,it,ir,jg)+dit*opl(ml,jt,ir,jg)) &
                      +dir*(djt*opl(ml,it,jr,jg)+dit*opl(ml,jt,jr,jg))
!	if in stretched layer get tgi from looking ahead?
  230           if(inwind) then
                  if(wmdot*vwind.ne.zero) &
                    dgw = wmdot/(vwind*r**2)
                  dens = dgw
                  pres = dens*tgi
                  ddfac = zero
                else
                  dens = pres/tgi
!	tgrav should change when H2 becomes H
                  hh2 = (tgi/tgrav)*(r**3/rstar)
!	dust density decrease exponential factor for a = 1 cm
                  if((avisc.gt.zero).and.(dens.gt.zero)) then
                    ddfac = ddfac+z*dz/(avisc*hh2*sqrt(hh2)*dens)
                  elseif((avisc.lt.zero).and.(dens.gt.zero)) then
                    avist = -avisc*min(ten,max(one,tgi/hun))**4
                    ddfac = ddfac+z*dz/(avist*hh2*sqrt(hh2)*dens)
                  else
                    ddfac = zero
                  endif
!	density change in half cell
                  gaus = sexp(-(z*dz)/(two*hh2))
!	pressure at center of cell
                  pres = pres*gaus
                  if(pres.lt.one) pres = zero
                  if((costh.gt.coswind).or. &
                    ((pres.lt.pram).and.(costh.ge.cosw))) then
                    inwind = .true.
                    cosw = max(cosw,costh-dz/r)
                    goto 230
                  endif
                  dens = pres/tgi
!	pressure at top of cell
                  pres = pres*gaus
                endif  !  inwind?
              endif  !  copy from inner grid
              tgz(iz,ip,ig) = tgi
              dgz(iz,ip,ig) = dens
              ddz(iz,ip,ig) = ddfac
!	leave taufac in for calculation of mdisk
              colup = colup+taufac*dens*dz
            enddo  !  iz
            xmdisk = xmdisk+fpi*rho*drho*colup
          enddo  !  ip
        enddo  !  igr

!	interpolate onto trg grid
!	could recalculate opl(ml) and condens here
        rphot = zero
        taucont = zero
        do igr = 0,ng
!	should always be in rho-z grid ig
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          drho = dr/four
          dz = dr/four
          lp = 4*lr
          kp = lp/2+1
          lz = 4*nr
          kz = lz/2+1
          do ir = kr,lr
            r = dr*(ir-half)
            colup = zero
            do it = 1,nt
              costh = dcosth*(nt-it+half)
              z = r*costh
              rho = sqrt(r**2-z**2)
              xiz = z/dz+half
              iz = xiz
              iz = max(1,min(lz-1,iz))
              jz = iz+1
              diz = xiz-iz
!	don''t extrapolate
              diz = max(zero,min(one,diz))
              djz = one-diz
              xip = rho/drho+half
              ip = xip
              ip = max(1,min(lp-1,ip))
              jp = ip+1
              dip = xip-ip
              dip = max(zero,min(one,dip))
              djp = one-dip
              dens = djp*(djz*dgz(iz,ip,ig)+diz*dgz(jz,ip,ig)) &
                    +dip*(djz*dgz(iz,jp,ig)+diz*dgz(jz,jp,ig))
              ddfac = djp*(djz*ddz(iz,ip,ig)+diz*ddz(jz,ip,ig)) &
                     +dip*(djz*ddz(iz,jp,ig)+diz*ddz(jz,jp,ig))
              dens = max(zero,dens)
              dg(it,ir,ig) = dens
              dd(it,ir,ig) = max(zero,ddfac)
              dm(it,ir,ig) = dens*(r/r0)**2
              if((ig.eq.ng9).and.(it.eq.1)) then
                hh2 = (tg(it,ir,ig)/tgrav)*(r**3/rstar)
                colup = dens*(hh2/z)
              endif
!	leave out opl(ml) since stretching is radial
              colup = colup+dens*r*dcosth/sqrt(one-costh**2)
              colz(it,ir,ig) = colup
            enddo  !  it
            dtau = (opa(iuv,nt,ir,ig)+sca(iuv,nt,ir,ig))*dr
            taucont = taucont+dtau
            if((rphot.eq.zero).and.(taucont.gt.one)) then
              dir = (taucont-one)/dtau
              rphot = dr*(ir+half-dir)
            endif
          enddo  !  ir
        enddo  !  ig
        call copyarr(dg)
        call copyarr(dm)
        call copyarr(dd)
        call copyarr(colz)

!	plot gas density
        if(doprint) print '(" disk mass:",1p,e10.2)', xmdisk*xmg/xmsun
        if(doplot) then
          print '(" plotting log_10(hydrogen density * r^2)")'
          monfile = ''
          if(nit.eq.(lit-1)) monfile = 'dens.con'
          call plotcon(dm,.true.,monfile)
          call waitasec(wait,2)
        endif
        if(ifmpi) call mpi_barrier(MPI_COMM_WORLD,ierr)

!	dust radiative transfer

        if(doprint) then
           print '(/" calculate dust radiative transfer  nit =",i3)', &
             nit
           write(16,'(" calculate dust radiative transfer  nit =",i3)') &
             nit
        endif

!	zero mean intensity arrays and spectrum
        call zerog(nu)
        do irot = 1,nrot
          do iu = 1,nu
            fw(iu,1,irot) = zero
          enddo
        enddo

!	loop over ray directions (grid orientations), farming out
      do irot0 = 1,nrot,mysize
        irot = irot0+myrank
        crot = dcrot*(nrot-irot+half)
        rot = acos(crot)

        call zerof(nu)

!	call outer grid routine, which loops over rays and calls inner grids
!	stepd is called from grids to do dust radiative transfer
!	and add contributions to mean intensity array
        call grid0(fwj,ny/2)

!	add contribution to spectrum from +y and save in fw
        fwmn = two*fwj(1)*dud(2)/fwmx
        do iu = 1,nu
          fwiu = two*fwj(iu)
          fw(iu,1,irot) = fwiu
          fwnorm = fwiu*dud(iu)/fwmx
          fwj(iu) = slog10(max(fwmn,fwnorm))
        enddo
        fwj(1) = slog10(fwmn)

!	broadcast calculated spectra and add contributions to mean intensity
        if(ifmpi) then
          do irank = 0,mysize-1
            call mpi_bcast(fw(1,1,irot0+irank),nu, &
            MPI_DOUBLE_PRECISION,irank,MPI_COMM_WORLD,ierr)
          enddo
          do igr = 0,ng
            call startgr(igr,ig,kr,lr,lt,dr,dcosth)
            do ir = kr,lr
              call mpi_allreduce(MPI_IN_PLACE,fn(1,ir,ig),nt, &
                MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
              call mpi_allreduce(MPI_IN_PLACE,ff(1,1,ir,ig),mu*nt, &
                MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
            enddo
          enddo
        endif
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            do it = 1,nt
              gn(it,ir,ig) = gn(it,ir,ig)+fn(it,ir,ig)
              do iu = 1,nu
                gg(iu,it,ir,ig) = gg(iu,it,ir,ig)+ff(iu,it,ir,ig)
              enddo
            enddo
          enddo
        enddo

!	plot images and spectrum
        if(doplot) then
          call movie(nc,1)
          call sleep(2)
          call movie(nc,2)
          call sleep(2)
          print '(" plotting log disk SED for irot =",i4)', irot
          jrot = (irot-1)/mysize+1
          call plotspec(fwj,mw,nu,1,jrot,zero,'')
        endif

!	end loop over grid orientations
      enddo  !  irot

!	fill in inner halves so I can smooth radially
        call copyarr(gn)
        call copyarrr(gg,mu,nu)
!	divide gg and qq by gn and average if appropriate
        call ggavg(nu,ntsm,nrsm)
        call copyarrr(gg,mu,nu)
!	plot gg before adding starlight
!        if(doplot) then
!          print '(" plotting disk light / stellar")'
!          call plotggfs(gg,ng9,iupl)
!          call waitasec(wait,2)
!        endif

!	radiative transfer of starlight (w/o sources)
!	note: not doing star for nt < irot < nnt
        if(doprint) print '(" starlight radiative transfer")'
!	first do nt rot steps for grids 1-ng
        call zerof(nu)
        dxstar = rstar/four
        idx = 4
        ix1 = nx9/2+1-idx
        ix2 = nx9/2+idx
!	set iy to a value unused by dust rad xfer
        iy = nx9/2+1
        do irot0 = 1,nt,mysize
          irot = irot0+myrank
          crot = dcost*(nt-irot+half)
          if(abs(crot).gt.one) goto 240
          rot = acos(crot)
          if((irot.ge.1).and.(irot.le.nt)) then
            do iu = 1,nu
              fw(iu,2,irot) = zero
            enddo
          endif
          sumbg = zero
          do ix = ix1,ix2
            x = dxstar*(ix-nx9/2-half)
            gwt = sqrt(one-(x/rstar)**2)
            sumbg = sumbg+bstar*gwt
!	start with stellar double blackbody spectrum
            do iu = 1,nu
              fwj(iu) = bus(iu)
            enddo
            do igr = 1,ng
              call startgr(igr,ig,kr,lr,lt,dr,dcosth)
              if(igr.eq.1) then
                lz = nz9
                dz = dz9
                kz = lz/2+1+nint(rstar/dz)
              else
                lz = nz
                kz = 3*lz/4+1
                twig = two**ig
                dz = dz0*twig
              endif
              do iz = kz,lz
                call stepd(fwj,ix,iy,iz,ig)
              enddo
            enddo  !  igr
!	sum outgoing stellar spectra in fw(iu,2,irot)
            if((irot.ge.1).and.(irot.le.nt)) then
              do iu = 1,nu
                fw(iu,2,irot) = fw(iu,2,irot)+bstar*gwt*fwj(iu)
              enddo
            endif
          enddo  !  ix
          if(ifmpi) then
!	can''t use allgather since arrays aren''t contiguous
!	but this caused computer to freeze if not all processors bcast
            do irank = 0,mysize-1
              irot = irot0+irank
              if((irot.ge.1).and.(irot.le.nt)) &
                call mpi_bcast(fw(1,2,irot),nu, &
                MPI_DOUBLE_PRECISION,irank,MPI_COMM_WORLD,ierr)
            enddo
          endif
  240   enddo  !  irot

!	do finer rot steps for grid 9
!	first zero ff in grid 9
        ig = ng9
        do ir = 1,nr9
          do it = 1,nt
            do iu = 1,nu
              ff(iu,it,ir,ig) = zero
            enddo
            fn(it,ir,ig) = zero
          enddo
        enddo
        igr = 0
        call startgr(igr,ig,kr,lr,lt,dr,dcosth)
        lz = nz9
        dz = dz9
        kz = lz/2+nint(rstar/dz+one)
        do irot0 = 1,nt,mysize
          irot = irot0+myrank
          crot = dcos9*(nt-irot+half)
          rot = acos(crot)
          do ix = ix1,ix2
!	start with stellar double blackbody spectrum
            do iu = 1,nu
              fwj(iu) = bus(iu)
            enddo
            do iz = kz,lz
              call stepd(fwj,ix,iy,iz,ig)
            enddo
          enddo  !  ix
        enddo  !  irot

!	add starlight to dust spectra and compare outgoing to stellar
        ffstar = zero
        fnstar = zero
        starsum = zero
        xtrot = real(nnt)/real(nrot)
        ntrot = nint(xtrot)
        do irot = 1,nrot
          it1 = nt-xtrot*(nrot-irot+one)+one+half
          it2 = nt-xtrot*(nrot-irot)+half
          starf = zero
          if(it1.ge.1) then
            do it = it1,it2
!	now include u < ud(1) by modifying dud(1)
!              starf = starf+fw(1,2,it)*ud(1)/(three*xtrot)
              do iu = 1,nu
                fw(iu,1,irot) = fw(iu,1,irot)+fw(iu,2,it)/xtrot
                starf = starf+fw(iu,2,it)*dud(iu)/xtrot
              enddo
            enddo
          else
            do iu = 1,nu
              fw(iu,1,irot) = fw(iu,1,irot)+astar*bus(iu)
              starf = starf+astar*bus(iu)*dud(iu)
            enddo
          endif
          sumfw = fw(1,1,irot)*ud(1)/four
          do iu = 1,nu
            sumfw = sumfw+fw(iu,1,irot)*dud(iu)
          enddo
          ffstar = ffstar+sumfw
          starsum = starsum+starf
          fnstar = fnstar+one
        enddo
        ffstar = ffstar/(fnstar*sumfs)
        starsum = starsum/(fnstar*sumfs)
        hcsum = hcsum/(xlstar*xlsun)
!        ffstar = ((ffstar+hcsum+one)/two-starsum)/(one-starsum)
        xlaccr = gmmdot/(two*rin*xlstar*xlsun)
        if(doprint) then
          print '(" outgoing flux/stellar:",f10.3)', ffstar
          print '(" just starlight:",f10.3)', starsum
          print '(" viscous heating, net gas heating power/stellar:", &
            1p,2e10.2)', xlaccr,hcsum
          write(16,'(" outgoing flux/stellar:",f10.3)') ffstar
          write(16,'(" just starlight:",f10.3)') starsum
        endif
        ffstar = (ffstar-starsum)/(one+hcsum-starsum)
        if((nit.lt.4).or.(nit.eq.lit-1)) then
          ffstar = one
        else
          ffstar = min(1.7,max(0.6,ffstar))
        endif
        if(doprint) print '(" correction factor:",f10.3)', one/ffstar

        if(doplot) then
!	plot SEDs including starlight
          print '(" plotting SEDs including starlight")'
          fwmn = fw(1,1,1)*dud(2)/fwmx
          do ipl = 1,8
            irot = (ipl*nrot)/8
            spec(1,ipl) = slog10(fwmn)
            do iu = 2,nu
              fwnorm = fw(iu,1,irot)*dud(iu)/fwmx
              spec(iu,ipl) = slog10(max(fwmn,fwnorm))
            enddo
          enddo
          if(nit.eq.lit) monfile = 'seds.spec'
          call plotspec(spec,mw,nu,9,1,zero,monfile)
          call waitasec(wait,2)
        endif
        if(doprint) then
          do ipl = 1,8
            irot = (ipl*nrot)/8
            xn13 =-slog10((ud(36)*fw(36,1,irot))/(ud(30)*fw(30,1,irot))) &
                  /slog10(ud(36)/ud(30))
            xn25 =-slog10((ud(30)*fw(30,1,irot))/(ud(26)*fw(26,1,irot))) &
                  /slog10(ud(30)/ud(26))
            xn31 =-slog10((ud(30)*fw(30,1,irot))/(ud(24)*fw(24,1,irot))) &
                  /slog10(ud(30)/ud(24))
!            print '(" SED indices:",3f8.3)', xn13,xn25,xn31
            if(nit.eq.lit) &
              write(16,'(" SED indices:",3f8.3)') xn13,xn25,xn31
          enddo
        endif

!	sum starlight mean intensities
        if(ifmpi) then
          do igr = 0,ng
            call startgr(igr,ig,kr,lr,lt,dr,dcosth)
            do ir = kr,lr
              call mpi_allreduce(MPI_IN_PLACE,fn(1,ir,ig),nt, &
                MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
              call mpi_allreduce(MPI_IN_PLACE,ff(1,1,ir,ig),mu*nt, &
                MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
            enddo
          enddo
        endif
!	add starlight to mean intensities
        nzero = 0
        sumfn = zero
        sumn = zero
        fnmin = one
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            r = dr*(ir-half)
            if(r.lt.rstar) then
              ostar = one
            elseif(r.lt.(hun*rstar)) then
              ostar = (one-sqrt(one-(rstar/r)**2))/two
            else
              ostar = quart*(rstar/r)**2
            endif
            do it = 1,nt
              if(smooth.and.(igr.eq.0)) then
!	3 point angular average in grid 9
!	is this a good idea?
                jt1 = it-1
                jt2 = it+1
                fni = zero
                do jt = jt1,jt2
                  if(jt.lt.1) then
                    kt = 1
                  elseif(jt.gt.nt) then
                    kt = 2*nt-jt+1
                  else
                    kt = jt
                  endif
                  fni = fni+fn(kt,ir,ig)
                enddo
              else
                fni = fn(it,ir,ig)
              endif
              if(fni.gt.fnmin) then
                sumfn = sumfn+fni
                sumn = sumn+one
                wtf = ostar/fni
                if(smooth.and.(igr.eq.0)) then
!	put ffstar correction into gg before adding starlight
                  do iu = 1,nu
                    gg(iu,it,ir,ig) = gg(iu,it,ir,ig)/ffstar
                  enddo
!	3 point angular smooth of ff in grid 9
                  jt1 = it-1
                  jt2 = it+1
                  do jt = jt1,jt2
                    if(jt.lt.1) then
                      kt = 1
                    elseif(jt.gt.nt) then
                      kt = 2*nt-jt+1
                    else
                      kt = jt
                    endif
                    do iu = 1,nu
                      gg(iu,it,ir,ig) = gg(iu,it,ir,ig) &
                        +wtf*ff(iu,kt,ir,ig)
                    enddo
                  enddo
                else
!	put ffstar correction into gg while adding starlight
                  do iu = 1,nu
                    gg(iu,it,ir,ig) = gg(iu,it,ir,ig)/ffstar &
                        +wtf*ff(iu,it,ir,ig)
                  enddo
                endif
              else
!	leave out ff?  could I use ff(ir+1)?
                do iu = 1,nu
                  gg(iu,it,ir,ig) = gg(iu,it,ir,ig)/ffstar
                enddo
                if((it.le.lt).and.(r.gt.rin)) nzero = nzero+1
              endif
            enddo  !  it
            do it = 1,nt
              if(nit.eq.(lit-1)) gguv(it,ir,ig) = gg(iufl,it,ir,ig)
              if((fn(it,ir,ig).gt.zero).and.(r.gt.two*rstar)) then
                do iu = 1,nu
                  ff(iu,it,ir,ig) = ostar*ff(iu,it,ir,ig)/fn(it,ir,ig)
                enddo
              else
                do iu = 1,nu
                  ff(iu,it,ir,ig) = zero
                enddo
              endif
            enddo
          enddo  !  ir
        enddo  !  igr
        if(doprint) &
          print '(" for starlight, mean rays/cell, zeros:",f8.2,i8)', &
            sumfn/sumn,nzero

!	plot mean intensity / stellar on xy grid
        call copyarrr(gg,mu,nu)
        if(doplot) then
          if(rim.lt.rgrid1) then
!            print '(" plotting starlight / unattentuated")'
!            call plotggfs(ff,ng9,iupl)
!            call waitasec(wait,2)
            print '(" plotting disk light + starlight / stellar")'
            call plotggfs(gg,ng9,iupl)
            call waitasec(wait,2)
            print '(" plotting disk light + starlight / stellar")'
            call plotggfs(gg,igpl,iupl)
            call waitasec(wait,2)
          endif
        endif

        if(doprint) &
          print '(" calculate gas and dust temperatures")'
!	first calculate temperatures in grid 9
!	loop over angles, farming out
        dtsum = zero
        dnsum = zero
        dtmax = zero
        igdtmx = 0
        irdtmx = 0
        itdtmx = 0
        hcsum = zero
        do it0 = 1,nt,mysize
          it = it0+myrank
          colin = zero
          taucont = zero
          igr = 0
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            r = dr*(ir-half)
            dens = dg(it,ir,ig)
!	remove exponential decay and inner rolloff from sigma
            sigma = two*colz(nt,ir,ig)/(one+(r0/r)**(qcol+pcol))
            if(r1.gt.zero) sigma = sigma*sexp(r/r1)
            sigma = max(sigma,1.0d+18)
            tgi = tg(it,ir,ig)
            tdi = td(it,ir,ig)
            pres = dens*tgi
            tsubli = tsubl*max(one,pres)**gsubl
            colup = colz(it,ir,ig)
            xwat = xwat0*sexp(-(tfrwat/tdi)**10)
            xco = xco0*sexp(-(tfrco/tdi)**10)
!	enhance CO where carbon is burned
            if((xburn.gt.one).and.(tburn.gt.one)) &
              xco = xco*(xburn+(one-xburn)*sexp(-(tgi/tburn)**10))
            if(taucont.le.zero) then
              xwat = xwat/hun
              xco = xco/hun
            endif
!	gastemp calculates temperatures of grain species and gas
!	from equations of heating and cooling
            call gastemp(it,ir,ig,iamax,tsubli,tgi,tdi,hcgas)
            dti = tgi-tg(it,ir,ig)
            dtsum = dtsum+dti**2
            dnsum = dnsum+one
            if(abs(dti).gt.abs(dtmax)) then
              dtmax = dti
              igdtmx = ig
              irdtmx = ir
              itdtmx = it
            endif
            dens = pres/tgi
            hcsum = hcsum+hcgas*dens*fpi*r**2*dr*dcosth
!	burn carbon grains; takes effect after next iteration, then permanent
!	but don''t burn carbon in inner disk
            if((tgi.gt.tburn).and.(r.gt.r0/two)) then
              do ia = 1,nac
                condens(ia,it,ir,ig) = zero
              enddo
            endif
            tg(it,ir,ig) = tgi
            td(it,ir,ig) = tdi
            dt(it,ir,ig) = dti
            dm(it,ir,ig) = tdi/tsubli
!            dm(it,ir,ig) = tdi/tgi
            if(pres.lt.one) dm(it,ir,ig) = one
!            dm(it,ir,ig) = sqrt(opa(iusil,it,ir,ig) &
!              *bw(ud(iusil),tdi)*r**3)
            colin = colin+opl(ml,it,ir,ig)*dens*dr
            colr(it,ir,ig) = colin
            taucont = taucont+(opa(iuv,it,ir,ig)+sca(iuv,it,ir,ig))*dr
            tauv(it,ir,ig) = taucont
          enddo  !  ir
!	broadcast after doing ir loop
          if(ifmpi) then
            do ir = kr,lr
              call mpi_allgather( &
                tg(it,ir,ig),1,MPI_DOUBLE_PRECISION, &
                tg(it0,ir,ig),1,MPI_DOUBLE_PRECISION, &
                MPI_COMM_WORLD,ierr)
              call mpi_allgather( &
                td(it,ir,ig),1,MPI_DOUBLE_PRECISION, &
                td(it0,ir,ig),1,MPI_DOUBLE_PRECISION, &
                MPI_COMM_WORLD,ierr)
              call mpi_allgather( &
                dt(it,ir,ig),1,MPI_DOUBLE_PRECISION, &
                dt(it0,ir,ig),1,MPI_DOUBLE_PRECISION, &
                MPI_COMM_WORLD,ierr)
              call mpi_allgather( &
                dm(it,ir,ig),1,MPI_DOUBLE_PRECISION, &
                dm(it0,ir,ig),1,MPI_DOUBLE_PRECISION, &
                MPI_COMM_WORLD,ierr)
              call mpi_allgather( &
                colr(it,ir,ig),1,MPI_DOUBLE_PRECISION, &
                colr(it0,ir,ig),1,MPI_DOUBLE_PRECISION, &
                MPI_COMM_WORLD,ierr)
              call mpi_allgather( &
                tauv(it,ir,ig),1,MPI_DOUBLE_PRECISION, &
                tauv(it0,ir,ig),1,MPI_DOUBLE_PRECISION, &
                MPI_COMM_WORLD,ierr)
              call mpi_allgather( &
                opa(1,it,ir,ig),mu,MPI_DOUBLE_PRECISION, &
                opa(1,it0,ir,ig),mu,MPI_DOUBLE_PRECISION, &
                MPI_COMM_WORLD,ierr)
              call mpi_allgather( &
                sca(1,it,ir,ig),mu,MPI_DOUBLE_PRECISION, &
                sca(1,it0,ir,ig),mu,MPI_DOUBLE_PRECISION, &
                MPI_COMM_WORLD,ierr)
              call mpi_allgather( &
                swd(1,it,ir,ig),mu,MPI_DOUBLE_PRECISION, &
                swd(1,it0,ir,ig),mu,MPI_DOUBLE_PRECISION, &
                MPI_COMM_WORLD,ierr)
              call mpi_allgather( &
                opl(1,it,ir,ig),ml,MPI_DOUBLE_PRECISION, &
                opl(1,it0,ir,ig),ml,MPI_DOUBLE_PRECISION, &
                MPI_COMM_WORLD,ierr)
              call mpi_allgather( &
                condens(1,it,ir,ig),ma,MPI_DOUBLE_PRECISION, &
                condens(1,it0,ir,ig),ma,MPI_DOUBLE_PRECISION, &
                MPI_COMM_WORLD,ierr)
            enddo  !  ir
          endif
        enddo  !  it0

!	calculate temperatures in grids 1-8
!	loop over angles, farming out
        do it0 = 1,nt,mysize
          it = it0+myrank
          costh = dcosth*(nt-it+half)
          lt9 = nt-nn9*(nt-it)
          kt9 = lt9-nn9+1
          colin = zero
          taucont = zero
          do igr = 1,ng
            call startgr(igr,ig,kr,lr,lt,dr,dcosth)
            if((igr.eq.1).and.(kt9.ge.1).and.(lt9.le.nt)) then
              do ir = kr,lr
                lr9 = ir*nn9
                kr9 = lr9-nn9+1
                sumtg = zero
                sumtd = zero
                sumcolr = zero
                sumtauv = zero
                do it9 = kt9,lt9
                  do ir9 = kr9,lr9
                    sumtg = sumtg+tg(it9,ir9,ng9)
                    sumtd = sumtd+td(it9,ir9,ng9)
                  enddo
                  sumcolr = sumcolr+colr(it9,lr9,ng9)
                  sumtauv = sumtauv+tauv(it9,lr9,ng9)
                enddo
                tg(it,ir,ig) = sumtg/xnn9
                td(it,ir,ig) = sumtd/xnn9
                colin = sumcolr/xn9
                colr(it,ir,ig) = colin
                taucont = sumtauv/xn9
                tauv(it,ir,ig) = taucont
              enddo
            else
              do ir = kr,lr
                r = dr*(ir-half)
                if(r.lt.rin) taucont = zero
                dens = dg(it,ir,ig)
                sigma = two*colz(nt,ir,ig)/(one+(r0/r)**(qcol+pcol))
                if(r1.gt.zero) sigma = sigma*sexp(r/r1)
                sigma = max(sigma,1.0d+18)
                tgi = tg(it,ir,ig)
                tdi = td(it,ir,ig)
                pres = dens*tgi
                tsubli = tsubl*max(one,pres)**gsubl
                colup = colz(it,ir,ig)
                xwat = xwat0*sexp(-(tfrwat/tdi)**10)
                xco = xco0*sexp(-(tfrco/tdi)**10)
!	enhance CO where carbon is burned
                if((xburn.gt.one).and.(tburn.gt.one)) &
                  xco = xco*(xburn+(one-xburn)*sexp(-(tgi/tburn)**10))
                if(taucont.le.zero) then
                  xwat = xwat/hun
                  xco = xco/hun
                endif
                call gastemp(it,ir,ig,iamax,tsubli,tgi,tdi,hcgas)
                dti = tgi-tg(it,ir,ig)
                if((ig.ne.1).and.(opa(iuv,it,ir,ig).gt.zero)) then
                  dtsum = dtsum+dti**2
                  dnsum = dnsum+one
                  if(abs(dti).gt.abs(dtmax)) then
                    dtmax = dti
                    igdtmx = ig
                    irdtmx = ir
                    itdtmx = it
                  endif
                endif
                dens = pres/tgi
                hcsum = hcsum+hcgas*dens*fpi*r**2*dr*dcosth
!                if(tgi.gt.tburn) then
                if((tgi.gt.tburn).and.(r.gt.r0/two)) then
                  do ia = 1,nac
                    condens(ia,it,ir,ig) = zero
                  enddo
                endif
                tg(it,ir,ig) = tgi
                td(it,ir,ig) = tdi
                dt(it,ir,ig) = dti
                dm(it,ir,ig) = tdi/tsubli
!                dm(it,ir,ig) = tdi/tgi
                if(pres.lt.one) dm(it,ir,ig) = one
!                dm(it,ir,ig) = sqrt(opa(iusil,it,ir,ig) &
!                  *bw(ud(iusil),tdi)*r**3)
                colin = colin+opl(ml,it,ir,ig)*dens*dr
                colr(it,ir,ig) = colin
                taucont = taucont &
                  +(opa(iuv,it,ir,ig)+sca(iuv,it,ir,ig))*dr
                tauv(it,ir,ig) = taucont
              enddo  !  ir
            endif  !  copy from grid 9?
          enddo  !  igr
          if(ifmpi) then
            do igr = 1,ng
              call startgr(igr,ig,kr,lr,lt,dr,dcosth)
              do ir = kr,lr
                call mpi_allgather( &
                  tg(it,ir,ig),1,MPI_DOUBLE_PRECISION, &
                  tg(it0,ir,ig),1,MPI_DOUBLE_PRECISION, &
                  MPI_COMM_WORLD,ierr)
                call mpi_allgather( &
                  td(it,ir,ig),1,MPI_DOUBLE_PRECISION, &
                  td(it0,ir,ig),1,MPI_DOUBLE_PRECISION, &
                  MPI_COMM_WORLD,ierr)
                call mpi_allgather( &
                  dt(it,ir,ig),1,MPI_DOUBLE_PRECISION, &
                  dt(it0,ir,ig),1,MPI_DOUBLE_PRECISION, &
                  MPI_COMM_WORLD,ierr)
                call mpi_allgather( &
                  dm(it,ir,ig),1,MPI_DOUBLE_PRECISION, &
                  dm(it0,ir,ig),1,MPI_DOUBLE_PRECISION, &
                  MPI_COMM_WORLD,ierr)
                call mpi_allgather( &
                  colr(it,ir,ig),1,MPI_DOUBLE_PRECISION, &
                  colr(it0,ir,ig),1,MPI_DOUBLE_PRECISION, &
                  MPI_COMM_WORLD,ierr)
                call mpi_allgather( &
                  tauv(it,ir,ig),1,MPI_DOUBLE_PRECISION, &
                  tauv(it0,ir,ig),1,MPI_DOUBLE_PRECISION, &
                  MPI_COMM_WORLD,ierr)
                call mpi_allgather( &
                  opa(1,it,ir,ig),mu,MPI_DOUBLE_PRECISION, &
                  opa(1,it0,ir,ig),mu,MPI_DOUBLE_PRECISION, &
                  MPI_COMM_WORLD,ierr)
                call mpi_allgather( &
                  sca(1,it,ir,ig),mu,MPI_DOUBLE_PRECISION, &
                  sca(1,it0,ir,ig),mu,MPI_DOUBLE_PRECISION, &
                  MPI_COMM_WORLD,ierr)
                call mpi_allgather( &
                  swd(1,it,ir,ig),mu,MPI_DOUBLE_PRECISION, &
                  swd(1,it0,ir,ig),mu,MPI_DOUBLE_PRECISION, &
                  MPI_COMM_WORLD,ierr)
                call mpi_allgather( &
                  opl(1,it,ir,ig),ml,MPI_DOUBLE_PRECISION, &
                  opl(1,it0,ir,ig),ml,MPI_DOUBLE_PRECISION, &
                  MPI_COMM_WORLD,ierr)
                call mpi_allgather( &
                  condens(1,it,ir,ig),ma,MPI_DOUBLE_PRECISION, &
                  condens(1,it0,ir,ig),ma,MPI_DOUBLE_PRECISION, &
                  MPI_COMM_WORLD,ierr)
              enddo  !  ir
            enddo  !  igr
          endif  !  mpi
        enddo  !  it
        if(ifmpi) then
          call mpi_allreduce(MPI_IN_PLACE,dtsum,1, &
            MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
          call mpi_allreduce(MPI_IN_PLACE,dnsum,1, &
            MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
          call mpi_allreduce(MPI_IN_PLACE,hcsum,1, &
            MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
        endif

        call copywta(tg,dg)
        call copywta(td,dg)
        call copyarr(dt)
        call copyarr(dm)
        call copytau(colr)
        call copytau(tauv)
        call copyarrr(opa,mu,nu)
        call copyarrr(sca,mu,nu)
        call copyarrr(swd,mu,nu)
        call copyarrr(ff,mu,nu)
        call copyarrr(gg,mu,nu)
        call copyarrr(opl,ml,ml)

        rmsdt = sqrt(dtsum/dnsum)
        if(doprint) then
          print '(" midplane sublimation, photosphere radii:", &
            2f10.4)', rim/au,rphot/au
          print '(" rms, max dT:",1p,2e10.2,3i6)', &
            rmsdt,dtmax,itdtmx,irdtmx,igdtmx
          write(16,'(" midplane sublimation, photosphere radii:", &
            2f10.4)') rim/au,rphot/au
          write(16,'(" rms, max dT:",1p,2e10.3,3i6)') &
            rmsdt,dtmax,itdtmx,irdtmx,igdtmx
        endif

        if(doplot) then
!	make plot of density and temperature
          do igr = 1,ng
            call startgr(igr,ig,kr,lr,lt,dr,dcosth)
            do ir = kr,lr
              r = dr*(ir-half)
              iw = ir+(ig-1)*nr/2
              if((ig.le.ng).and.(iw.gt.0).and.(iw.le.mw)) then
                if(r.lt.rstar) r = rstar
                dens = dg(nt,ir,ig)
                colup = colz(nt,ir,ig)*xmg
                tgi = tg(nt,ir,ig)
                tdi = td(nt,ir,ig)
                tsubli = tsubl*max(one,dens*tgi)**gsubl
                do it = 1,nt
                  if(tauv(it,ir,ig).lt.one) itrim = it
                enddo
                tgrim = tg(itrim,ir,ig)
                tdrim = td(itrim,ir,ig)
                smin = 1.0d-02
                smax = 1.0d+06
                spec(iw,1) = slog10(max(smin,min(smax,colup)))
                spec(iw,2) = slog10(max(smin,min(smax,tgi)))
                spec(iw,3) = slog10(max(smin,min(smax,tdi)))
                spec(iw,4) = slog10(max(smin,min(smax,tgrim)))
                spec(iw,5) = slog10(max(smin,min(smax,tdrim)))
                spec(iw,6) = spec(iw,16)
              endif
            enddo  !  ir
          enddo  !  igr
          nwspec = min(mw,iw+16)
          do iw = nwspec-15,nwspec
            spec(iw,1) = -1.5
            spec(iw,2) = -0.5
            spec(iw,3) = 0.5
            spec(iw,4) = 1.5
            spec(iw,5) = 2.5
            spec(iw,6) = 3.5
          enddo

          monfile = ''
          if(nit.eq.lit) monfile = 'midplane'
          print '(" plotting log(Ng, mid Tg, Td, surf Tg, Td, T_C&G)")'
          call plotspec(spec,mw,nwspec,6,1,zero,monfile)
          call sleep(4)
!	plot tauv
          print '(" plotting log_10(tauv)")'
          call plotcon(tauv,.true.,'')
          call waitasec(wait,2)
!	plot temperature
          print '(" plotting log_10(gas temperature)")'
          if(nit.eq.lit) monfile = 'temp.con'
          call plotcon(tg,.true.,monfile)
          call waitasec(wait,2)
!	plot delta temperature
          print '(" plotting temperature change")'
          call plotcon(dt,.false.,'')
          call waitasec(wait,2)
!	plot dust temperature / sublimation
!          print '(" plotting dust temperature / sublimation")'
!          print '(" plotting silicate emission * r^3")'
!          call plotcon(dm,.false.,'')
!          call waitasec(wait,2)
        endif
        if(ifmpi) call mpi_barrier(MPI_COMM_WORLD,ierr)

!	loop back to hydro equil calculation
        if(nit.lt.lit) goto 200

!	save results of dust calculation
        if(doprint) then
          call savemodel
          open(unit=10,file='seds.mon',access='sequential')
          do iu = 1,nu
            write(10,'(i6,9f8.4)') iu,(spec(iu,ipl),ipl=1,9)
          enddo
          close(unit=10)
        endif

!	start radiative transfer in molecular lines
!	start here if reusing dust disk model
  300   if(ifmpi) call mpi_barrier(MPI_COMM_WORLD,ierr)
        moltype = abs(moltype)
        if(moltype.eq.0) then
          if(ifmpi) call mpi_finalize(ierr)
          stop
        endif

!	calculate colin and tauv.  probably not needed unless reuse = true
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            r = dr*(ir-half)
            do it = nt,1,-1
              if(ir.gt.kr) then
                colr(it,ir,ig) = colr(it,ir-1,ig) &
                  +opl(ml,it,ir-1,ig)*dg(it,ir-1,ig)*dr
                tauv(it,ir,ig) = tauv(it,ir-1,ig) &
                  +opa(iuv,it,ir-1,ig)*dr
              elseif(igr.gt.1) then
                colr(it,ir,ig) = colr(it,lr,ig-1) &
                  +opl(ml,it,lr,ig-1)*dg(it,lr,ig-1)*dr/two
                tauv(it,ir,ig) = tauv(it,lr,ig-1) &
                  +opa(iuv,it,lr,ig-1)*dr/two
              else
                colr(it,ir,ig) = zero
                tauv(it,ir,ig) = zero
              endif
            enddo
          enddo
        enddo
        call copytau(colr)
        call copytau(tauv)

!	set up velocity arrays
!	could add turbulence and radial pressure gradient influence
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          if(igr.le.1) itwind = nt
          do ir = kr,lr
            r = dr*(ir-half)
            if(wmdot*vwind.ne.zero) then
              dgw = 1.001d+00*wmdot/(vwind*r**2)
            else
              dgw = zero
            endif
            inwind = .false.
            do it = nt,1,-1
              costh = dcosth*(nt-it+half)
              if((it.le.itwind).and.(dg(it,ir,ig).lt.dgw)) then
                inwind = .true.
                itwind = it
              endif
!	wr and wth may be hardwired = 0
              if(r.lt.rstar) then
                wr(it,ir,ig) = zero
                wth(it,ir,ig) = zero
                wph(it,ir,ig) = vorb*sqrt(au/rstar)*sqrt(one-costh**2)
              elseif(inwind) then
                wr(it,ir,ig) = vwind
                wth(it,ir,ig) = zero
                wph(it,ir,ig) = zero
              else
                wr(it,ir,ig) = vterm*sqrt(one-rstar/r)
                wth(it,ir,ig) = zero
                wph(it,ir,ig) = vorb*sqrt(au/r)*sqrt(one-costh**2)
              endif
            enddo  !  it
          enddo  !  ir
        enddo  !  igr
        call copyarr(wr)
        call copyarr(wth)
        call copyarr(wph)

        nx = lx
        nrot = lrot
        dcrot = one/nrot
        ny = nx
        if(bipolar) then
          ly = lx/2
        else
          ly = lx
        endif
        nx2 = nx/2
        ny2 = ny/2
        nz2 = nz/2
        nx9 = nn9*nx
!	keep nz from dust model
        nz9 = nn9*nz
        xnx = nx
        xnz = nz
        dx0 = (two*nr*dr0)/xnx
        dx9 = two*dx0/xn9
        dz0 = xnx*dx0/xnz
        dz9 = two*dz0/xn9
        ll = nl
        lw = nw
!	first J,K
        jfirst = 50./sqrt(be)
        jlast = 1
        ju = jfirst
        if(moltype.eq.3) then
          kk = ju
        else
          kk = 0
        endif
        jr = ju-1

!	calculate molecular density etc.
        rphoto = zero
        rhi = zero
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            r = dr*(ir-half)
            do it = nt,1,-1
              tgi = max(ten,tg(it,ir,ig))
              tdi = max(ten,td(it,ir,ig))
              dens = dg(it,ir,ig)
              colin = colr(it,ir,ig)
!	assume molecules are destroyed by UV inside coldiss
!	and by heat above tdiss
              if(((colin.lt.coldiss).and.(tauv(it,ir,ig).lt.fifth)) &
                .or.(tgi.gt.tdiss)) then
                dmol = zero
                if(it.eq.nt) rphoto = r
!	assume molecules can''t form without dust
              elseif(tauv(it,ir,ig).eq.zero) then
                dmol = zero
                if(it.eq.nt) rphoto = r
              else
                dmol = dens*xmol*sexp(-(tfrmol/tdi)**10)
              endif
!	enhance CO where carbon is burned
              if((moltype.eq.1).and.(xburn.gt.one).and.(tburn.gt.one)) &
                dmol = dmol*(xburn+(one-xburn)*sexp(-(tgi/tburn)**10))
!	put optical depth decrease factor into dmol
              dmol = opl(ml,it,ir,ig)*dmol
              bkt = bk*tgi
              hcbkt = hcb/bkt
              qrot = bkt/hcb
!	ignore qvib?
!	modify qrot for symmetric tops
              if(hcc.gt.zero) qrot = (4./3.)*qrot*sqrt(bkt/hcc)
              if((ig.eq.ng).and.(ir.eq.nr)) dmol = zero
              qrot = bk*tgi/hcb
              dmq = dmol/qrot
!	save dmol/qrot in dm
              dm(it,ir,ig) = dmq
            enddo  !  it
          enddo  !  ir
        enddo  !  igr
        call copyarr(dm)
        if(doprint) print '(" molecules are destroyed inside of ", &
          f8.3," AU")',  rphoto/au

  400   if(doprint) then
          print '(/" Starting calculation for Ju,K =",2i3)', ju,kk
          write(16,'(/" Starting calculation for Ju,K =",2i3)') ju,kk
        endif
        nit = 1
!	initialize variables dependent on ju
        call initju(ju,kk)
        jr = ju-1
        gju = gnucl*(two*xju+one)
        gjr = gnucl*(two*xjr+one)
        eju = hcb*xju*(xju+one)
        ejr = hcb*xjr*(xjr+one)
        if(hcc.gt.zero) then
          ekk = (hcc-hcb)*xk*xk
        else
          ekk = zero
        endif
        if(lte) then
          dcru = zero
          dcrr = zero
        else
          dcru = (xju**4/gju)*arot/grot
          dcrr = (xjr**4/gjr)*arot/grot
        endif
        wj = two*be*xju
        wjr = two*be*xjr
        iuj = 16.*(slog10(wj)-one)
        iujr = 16.*(slog10(wjr)-one)
        iuj = max(2,min(nu,iuj))
        iujr = max(2,min(nu,iujr))
!	assume a flat stellar spectrum near each line
        fstar(1) = bw(wj,tstar)

        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            do it = 1,lt
              tgi = max(ten,tg(it,ir,ig))
              bkt = bk*tgi
              dmq = dm(it,ir,ig)
              dens = max(one,dg(it,ir,ig))
!	could try to include radiative excitation if I know gg
              if(lte) then
                dl(1,it,ir,ig) = dmq*gjr*sexp(-(ejr+ekk)/bkt)
                dl(2,it,ir,ig) = dmq*gju*sexp(-(eju+ekk)/bkt)
              elseif(reuse.and.(nit.eq.1)) then
                dcrit = dcrr/sqrt(tgi)
                dl(1,it,ir,ig) = dmq*gjr*sexp(-(ejr+ekk)/bkt) &
                                 *dens/(dens+dcrit)
                dcrit = dcru/sqrt(tgi)
                dl(2,it,ir,ig) = dmq*gju*sexp(-(eju+ekk)/bkt) &
                                 *dens/(dens+dcrit)
              elseif(nit.eq.1) then
!	try gg(iuj) if not reusing dust calc.
                ggiuj = gg(iujr,it,ir,ig)
                expbb = (ggiuj/(two*hc2*wjr**3+ggiuj))**(xju/two)
                if((expbb.lt.zilch).or.badnum(expbb)) expbb = zero
                dcrit = dcrr/sqrt(tgi)
                dl(1,it,ir,ig) = dmq*gjr*(dens*sexp(-(ejr+ekk)/bkt) &
                                 +dcrit*expbb)/(dens+dcrit)
                ggiuj = gg(iuj,it,ir,ig)
                expbb = expbb*(ggiuj/(two*hc2*wj**3+ggiuj))
                dcrit = dcru/sqrt(tgi)
                dl(2,it,ir,ig) = dmq*gju*(dens*sexp(-(eju+ekk)/bkt) &
                                 +dcrit*expbb)/(dens+dcrit)
              else
!	gg(1) has the mean intensity of the current line after the first nit
!	although that''s not really what I want
                ggj = gg(1,it,ir,ig)
!	I''d rather use the previous line here
                expbb = (ggj/(two*hc2*wj**3+ggj))
                dcrit = dcrr/sqrt(tgi)
                dl(1,it,ir,ig) = dmq*gjr*(dens*sexp(-(ejr+ekk)/bkt) &
                                 +dcrit*expbb)/(dens+dcrit)
                dcrit = dcru/sqrt(tgi)
                dl(2,it,ir,ig) = dmq*gju*(dens*sexp(-(eju+ekk)/bkt) &
                                 +dcrit*expbb)/(dens+dcrit)
              endif
            enddo  !  it
          enddo  !  ir
        enddo  !  igr
        call copyarrr(dl,ml,2)

!	radiative transfer along rays
        call zerog(1)

!	loop over ray directions (grid orientations)
      do irot0 = 1,nrot,mysize
        irot = irot0+myrank
        crot = dcrot*(nrot-irot+half)
        rot = acos(crot)
        call zerof(1)

!	call outer grid routine, which loops over rays and calls inner grids

        call grid0(fwj,ly)

        if(doprint) print &
          '(" grid orientation",i4," of",i4,"  rot =",f8.3)', &
          irot0+mysize/2,nrot,acos(dcrot*(nrot-irot0-mysize/2+half))
        if(doplot) then
          call movie(nc,1)
          call sleep(2)
        endif

!	add contribution to spectrum from +y (with opposite Doppler shifts)
!	and save and broadcast spectrum
        do iw = 1,nw
          if(bipolar.and.(iw.le.nw2)) then
            iw1 = nw+1-iw
            fwj(iw) = fwj(iw)+fwj(iw1)
            fwj(iw1) = fwj(iw)
          endif
          fw(iw,1,irot) = fwj(iw)
        enddo
        if(ifmpi) then
          do irank = 0,mysize-1
            call mpi_bcast(fw(1,1,irot0+irank),nw, &
            MPI_DOUBLE_PRECISION,irank,MPI_COMM_WORLD,ierr)
          enddo
        endif
        if(doplot) call plotspec(fwj,mw,nw,1,1,-dv/vkms,'')

!	add contributions to mean intensity (for possible non-LTE calc)
!	currently ff is not calculated
!        if(ifmpi) then
!          do igr = 0,ng
!            call startgr(igr,ig,kr,lr,lt,dr,dcosth)
!            do ir = kr,lr
!              call mpi_allreduce(MPI_IN_PLACE,ff(1,1,ir,ig),mu*nt, &
!                MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
!              call mpi_allreduce(MPI_IN_PLACE,fn(1,ir,ig),nt, &
!                MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
!            enddo
!          enddo
!        endif
!        do igr = 0,ng
!          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
!          do ir = kr,lr
!            do it = 1,nt
!              gn(it,ir,ig) = gn(it,ir,ig)+fn(it,ir,ig)
!              gg(ju,it,ir,ig) = gg(ju,it,ir,ig)+ff(1,it,ir,ig)
!            enddo
!          enddo
!        enddo

!	end loop over grid orientations
      enddo  !  irot

!	fill in inner halves so I can smooth radially
!        call copyarr(gn)
!        call copyarrr(gg,mu,1)
!	divide gg by gn and average if appropriate
!        call ggavg(1,ntsm,nrsm)

!	store final spectra

!	this would set the stellar continuum to 1
!        fnorm = fstar(1)*astar
!	this is supposed to give me Jy at 100 pc
        fnorm = 1.0d-23*clt*(3.086d+20/(dx0*two**ng))**2

        if(doplot) then
          mpl = (npl+7)/8
          w0 = (nw+1)/2.
          write(11,'("J =",3i3)') ju,jr,kk
          do jpl = 1,mpl
            lpl = min(npl,8*jpl)
            kpl = max(1,lpl-7)
            do ipl = kpl,lpl
              irot = (ipl*nrot)/npl
              call convo(fw(1,1,irot),spec(1,ipl), fnorm,dv/vkms,wj,dsm,nsm)
            enddo
            if(spec(1,1).lt.(5.0d+02)) then
              do iw = 1,nw
                vi = -dv*(iw-w0)/vkms
                write(11,'(i4,f8.2,8f8.3)') iw,vi,(spec(iw,ipl),ipl=kpl,lpl)
              enddo
            elseif(spec(1,1).lt.(5.0d+03)) then
              do iw = 1,nw
                vi = -dv*(iw-w0)/vkms
                write(11,'(i4,f8.2,8f8.2)') iw,vi,(spec(iw,ipl),ipl=kpl,lpl)
              enddo
            elseif(spec(1,1).lt.(5.0d+04)) then
              do iw = 1,nw
                vi = -dv*(iw-w0)/vkms
                write(11,'(i4,f8.2,8f8.1)') iw,vi,(spec(iw,ipl),ipl=kpl,lpl)
              enddo
            else
              do iw = 1,nw
                vi = -dv*(iw-w0)/vkms
                write(11,'(i4,f8.2,8f8.0)') iw,vi,(spec(iw,ipl),ipl=kpl,lpl)
              enddo
            endif
            call plotspec(spec(1,kpl),mw,nw,8,jpl,-dv/vkms,'')
            call waitasec(wait,2)
          enddo
        endif
        if(doprint) then
          write(16,'(" J =",3i3," spectra at 100 pc")') ju,jr,kk
          write(16, &
            '("      Fc(Jy) Fl(Jy-km/s)    Fl(W/m2)    EW(km/s)")')
          print '(" J =",3i3," spectra at 100 pc")', ju,jr,kk
          print '("      Fc(Jy) Fl(Jy-km/s)    Fl(W/m2)    EW(km/s)")'
          do ipl = 1,npl
            irot = (ipl*nrot)/npl
            call convo(fw(1,1,irot),spec(1,ipl), fnorm,dv/vkms,wj,dsm,nsm)
          enddo
        endif
        if(ifmpi) call mpi_barrier(MPI_COMM_WORLD,ierr)
  
!	could do isotopic and hot-band spectra
  
        if(kk.gt.1) then
          kk = kk-2
          goto 400
        endif
!        if(ju.le.8) then
!          ju = ju-2
!        elseif(ju.le.16) then
!          ju = ju-2
!        elseif(ju.le.32) then
!          ju = ju-4
!        else
!          ju = ju-8
!        endif
        ju = ju-2
        if(moltype.eq.3) kk = ju

        if(ju.ge.jlast) goto 400

        if(doplot) close(unit=11)
        if(doprint) close(unit=16)
        if(ifmpi) call mpi_finalize(ierr)

      end



      subroutine startgr(igr,ig,kr,lr,lt,dr,dcosth)

        implicit real*8 (a-h,o-z)
        parameter (two = 2.0d+00)
        save dr0,dr9,dcost,dcos9
        save nr,nr3,nr9,ng9,nt,nt8

        if((igr.eq.0).or.(igr.eq.ng9)) then
          ig = ng9
          dr = dr9
          kr = 1
          lr = nr9
          dcosth = dcos9
        else
          ig = igr
          dr = dr0*two**ig
          kr = nr3
          lr = nr
          dcosth = dcost
        endif
        if(ig.eq.1) then
          kr = 1
          lt = nt8
        else
          lt = nt
        endif
        return

      entry initgr(xr0,xr9,xcost,xcos9,ir,ir9,ig9,it,it8)

        dr0 = xr0
        dr9 = xr9
        dcost = xcost
        dcos9 = xcos9
        nr = ir
        nr9 = ir9
        ng9 = ig9
        nr3 = nr/2+1
        nt = it
        nt8 = it8
        return

      end



      subroutine dustinit(xc,xsi,ac,asi,tuv,xuv,emiss,tsubli)

        implicit real*8 (a-h,o-z)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mg = 20)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /angls/ dcost,dcos9,rot,irot
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu), &
                phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg), &
                opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg), &
                rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml), &
                wl(ml),dl(ml,mt,mr,mg),sl(ml,mt,mr,mg),opl(ml,mt,mr,mg)
        common /taus/ taucont,tauline,taucent,iuv,iusil,iufl,iwterm
        dimension spec(mu,ma),emis(mu)
        logical doprint,badnum
        parameter (pi = 3.14159265d+00)
        parameter (tpi = 6.28318530d+00)
        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (quart = 0.25d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (four = 4.0d+00)
        parameter (ten = 1.0d+01)
        parameter (hun = 1.0d+02)
        parameter (xmg = 4.0d-24)

        nr2 = nr/2
        doprint = (myrank.eq.0)
!	a = 0.3 nm - 100 um (may be changed)
!	na = 24; 12 each for carbon and silicates
        la = na/2
        arat = ten**half
        a0 = 1.0d-07
!	roll off distributions at ac and asi
!	grain distribution power laws, dn/dlna ~ a^-pow
        powc = 2.5d+00
        pows = 2.0d+00
!	assume AC grain density = 2 g/cc or 1e23 C atoms/cc
!	so grain vol/C = 1e-23, or vol/H = 1e-23*XC
!	solar XC = 3.6e-04, gas particles / H atom = 0.6
        xc = 3.6d-04*xc
        cvolg = (1.0d-23)*xc/(0.6d+00)
!	assume Sil grain density = 3.6 g/cc, with 172 AMU/Si,
!	or 1.25e22 Si/cc, or 8e-23 cc/Si
        xsi = 3.5d-05*xsi
        svolg = (8.0d-23)*xsi/(0.6d+00)
!	u = 11.5 - 10^5 cm-1
        u0 = ten
        dlogu = 4.0d+00/nu
        urat = ten**dlogu
        durat = (urat-one)/sqrt(urat)
!	w = 328 - 3162 cm-1
        w0 = ten**2.5
        dlogw = one/nu
        wrat = ten**dlogw
        dwrat = (wrat-one)/sqrt(wrat)
        wv = 2.0d+04
        iuv = slog10(wv/u0)/dlogu+half
        wfl = 5.0d+04
        iufl = slog10(wfl/u0)/dlogu+half
        wsil = 1.0d+03
        iusil = slog10(wsil/u0)/dlogu+half
        iwsil = slog10(wsil/w0)/dlogw+half
        if(iwsil.ne.iusil) &
          print '("Warning: iwsil.ne.iusil",2i6)',iwsil,iusil
        sumv = zero
        if(na.ne.24) print '("Warning: na =",i4)', na
        volc = zero
        vols = zero
        a = a0
        do ia = 1,la
          iac = ia
          ias = ia+la
          a = a*arat
          ad(iac) = a
          ad(ias) = a
!	approximation to Weingartner & Draine (2001)
!	note: popd = dn/dlna = a*dn/da since a is spaced logarithmically
          pops = one/((a/asi)**pows+(a/asi)**10)
          popc = one/((a/ac)**powc+(a/ac)**10)
          popd(iac) = popc
          popd(ias) = pops
          vola = (4.*pi/3.)*a**3
          volc = volc+popc*vola
          vols = vols+pops*vola
        enddo
        cnorm = cvolg/volc
        snorm = svolg/vols
        do ia = 1,la
          iac = ia
          ias = ia+la
          popd(iac) = cnorm*popd(iac)
          popd(ias) = snorm*popd(ias)
        enddo
        u = u0
        do iu = 1,nu
          u = u*urat
          ud(iu) = u
          dud(iu) = u*durat
          bus(iu) = bw(u,tstar)+xuv*(tstar/tuv)**4*bw(u,tuv)
        enddo
!	modify dud(1) so integral over ud includes ud < ud(1)
        dud(1) = ud(1)/four-dud(1)/two
!	initialize phab, phsc, reading from ac.out and as.out
        open(unit=10,file='ac.out',access='sequential')
        read(10,'(1x)')
        do ia = 1,la
          a = ad(ia)
!	skip every other a
          do iu = nu,1,-1
            read(10,'(///)')
          enddo
          do iu = nu,1,-1
!	read every 4th line
            read(10,'(6e12.4///)') amie,wlen,qext,qabs,qscat,cosscat
            if(doprint.and.((abs(1.0e+04-amie/a).gt.one) &
              .or.(abs(1.0e+04-ud(iu)*wlen).gt.one))) &
              print '(2i6,1p,4e10.2)', ia,iu,a,amie,ud(iu),wlen
            phab(iu,ia) = qabs
!	decrease qscat to account for lowered effect of forward scattering
            if((cosscat.gt.zero).and.(cosscat.lt.one)) &
              qscat = (one-cosscat)*qscat
            phsc(iu,ia) = qscat
          enddo
        enddo
        close(unit=10)
        open(unit=10,file='as.out',access='sequential')
        read(10,'(1x)')
        do ia = la+1,na
          a = ad(ia)
          do iu = nu,1,-1
            read(10,'(///)')
          enddo
          do iu = nu,1,-1
            read(10,'(6e12.4///)') amie,wlen,qext,qabs,qscat,cosscat
            if(doprint.and.((abs(1.0e+04-amie/a).gt.one) &
              .or.(abs(1.0e+04-ud(iu)*wlen).gt.one))) &
              print '(2i6,1p,4e10.2)', ia,iu,a,amie,ud(iu),wlen
            phab(iu,ia) = qabs
            if((cosscat.gt.zero).and.(cosscat.lt.one)) &
              qscat = (one-cosscat)*qscat
            phsc(iu,ia) = qscat
          enddo
        enddo
        close(unit=10)
!	finer grid for last iteration
!	neglect scattering.  it should be small beyond 3um?
        w = w0
        do iu = 1,nu
          w = w*wrat
          wd(iu) = w
        enddo
        open(unit=10,file='ac.out',access='sequential')
        read(10,'(1x)')
        do ia = 1,la
          a = ad(ia)
          do iu = nu,1,-1
            read(10,'(///)')
          enddo
!	skip to 31.6um
          do iu = 1,96
            read(10,'(1x)')
          enddo
          do iu = nu,1,-1
            read(10,'(5e12.4)') amie,wlen,qext,qabs,qscat
            if(doprint.and.((abs(1.0e+04-amie/a).gt.one) &
              .or.(abs(1.0e+04-wd(iu)*wlen).gt.one))) &
              print '(2i6,1p,4e10.2)', ia,iu,a,amie,wd(iu),wlen
            phid(iu,ia) = qabs
          enddo
          do iu = 1,96
            read(10,'(1x)')
          enddo
        enddo
        close(unit=10)
        open(unit=10,file='as.out',access='sequential')
        read(10,'(1x)')
        do ia = la+1,na
          a = ad(ia)
          do iu = nu,1,-1
            read(10,'(///)')
          enddo
          do iu = 1,96
            read(10,'(1x)')
          enddo
          do iu = nu,1,-1
            read(10,'(5e12.4)') amie,wlen,qext,qabs,qscat
            if(doprint.and.((abs(1.0e+04-amie/a).gt.one) &
              .or.(abs(1.0e+04-wd(iu)*wlen).gt.one))) &
              print '(2i6,1p,4e10.2)', ia,iu,a,amie,wd(iu),wlen
            phid(iu,ia) = qabs
          enddo
          do iu = 1,96
            read(10,'(1x)')
          enddo
        enddo
        close(unit=10)

!	calculate emis = 1-albedo averaged over stellar spectrum
!	but note that sqrt(avg) != avg(sqrt)
        emissum = zero
        emissil = zero
        sumbus = zero
        do iu = 1,nu
          busdud = bus(iu)*dud(iu)
          opac = zero
          scat = zero
          opacsil = zero
          scatsil = zero
          do ia = 1,na
            area = pi*ad(ia)**2
            popsig = popd(ia)*area*phab(iu,ia)
            if(mod(ia,2).eq.0) &
              spec(iu,ia/2) = max(-1.0d+01,(slog10(popsig)+2.4d+01))
            opac = opac+popsig
            if(ia.gt.la) opacsil = opacsil+popsig
            popsig = popd(ia)*area*phsc(iu,ia)
            scat = scat+popsig
            if(ia.gt.la) scatsil = scatsil+popsig
          enddo
          spec(iu,1) = max(-1.0d+01,(slog10(opac)+2.4d+01))
          emis(iu) = opac/(opac+scat)
          emissum = emissum+busdud*opac/(opac+scat)
          emissil = emissil+busdud*opacsil/(opacsil+scatsil)
          sumbus = sumbus+busdud
        enddo
        emiss = emissum/sumbus
        emissil = emissil/sumbus
        if(doprint) print '(" dust albedos and reflectivities:",4f8.3)', &
          one-emiss,one-emissil,one-sqrt(emiss),one-sqrt(emissil)
        treem = tsubli/two**quart
        reemiss = half*sqrt(emiss)*(treem/tstar)**4
        do iu = 1,nu
          rebus(iu) = bus(iu) &
            +sqrt(emis(iu))*(reemiss*bw(ud(iu),treem)-bus(iu))
        enddo

!	plot opacity and populations
        if(doprint) then
          tauv = zero
          tausil = zero
          taufir = zero
          wfir = 1.0d+02
          iufir = slog10(wfir/u0)/dlogu+half
          do ia = 1,na
            tauv = tauv+popd(ia)*pi*ad(ia)**2*phab(iuv,ia)
            tausil = tausil+popd(ia)*pi*ad(ia)**2*phab(iusil,ia)
            taufir = taufir+popd(ia)*pi*ad(ia)**2*phab(iufir,ia)
          enddo
          print '(" Tauv/N =",1p,4e10.2)', &
            tauv,tauv/xmg,tausil/xmg,taufir/xmg
          call plotspec(spec,mu,nu,na/2,1,zero,'')
          call sleep(4)
          do iu = 1,nu
            do ia = 2,na,2
              spec(iu,ia/2) = max(-5.,slog10(phab(iu,ia)))
            enddo
          enddo
          call plotspec(spec,mu,nu,na/2,1,zero,'')
        endif

        return
      end



      subroutine stepd(fwi,ix,iy,iz,igr)

        implicit real*8 (a-h,o-z)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mg = 20)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        parameter (mw = 1024)

        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg), &
                ff(mu,mt,mr,mg),fn(mt,mr,mg)
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /angls/ dcost,dcos9,rot,irot
        common /xj/ xju,xjp,xjq,xjr,xk
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu), &
                phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg), &
                opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg), &
                rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml), &
                wl(ml),dl(ml,mt,mr,mg),sl(ml,mt,mr,mg),opl(ml,mt,mr,mg)
        common /vels/ wr(mt,mr,mg),wth(mt,mr,mg),wph(mt,mr,mg)
        common /taus/ taucont,tauline,taucent,iuv,iusil,iufl,iwterm
        common /exptab/ expi(10239),too
        dimension fwi(mw),fwk(mw),ffi(ml)
        logical good,doprint,badnum
        logical starlt,newcell

        parameter (pi = 3.14159265d+00)
        parameter (sqrtpi = 1.772454d+00)
        parameter (sqrt2 = 1.414214d+00)
        parameter (alog2 = 0.6931472d+00)
        parameter (small = 1.0d-03)
        parameter (tiny = 1.0d-08)
        parameter (zilch = 1.0d-200)
        parameter (hun = 1.0d+02)
        parameter (big = 5.0d+03)
        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (four = 4.0d+00)
        parameter (eighth = 0.125d+00)
        parameter (au = 1.496d+13)

        save kt,kr,kg

        starlt = (iy.gt.nx9/2)
        jg = abs(igr)
        if(jg.le.1) jg = ng9
        if(starlt) then
          lx = nx9
          dx = rstar/four
          if(jg.eq.ng9) then
            jg = ng9
            lz = nz9
            dz = dz9
            rmx = nr9*dr9
          else
            twig = two**jg
            lz = nz
            dz = dz0*twig
            rmx = nr*dr0*twig
          endif
        elseif(jg.eq.ng9) then
          jg = ng9
          lx = nx9
          dx = dx9
          lz = nz9
          dz = dz9
          rmx = nr9*dr9
        else
          twig = two**jg
          lx = nx
          dx = dx0*twig
          lz = nz
          dz = dz0*twig
          rmx = nr*dr0*twig
        endif
        doprint = (myrank.eq.0).and.(ix.eq.1).and.(iy.eq.1) &
          .and.(jg.eq.ng9).and.(abs(iz-lz/2).lt.lz/8)
        ly = lx/2
!	finer xy steps in middle of grid
        if(igr.lt.0) dx = dx/two
        dx2 = dx/two
        dz2 = dz/two
        ds = dz
        x = dx2*(2*ix-lx-1)
        y = dx2*(2*iy-lx-1)
        z = dz2*(2*iz-lz-1)
        r = sqrt(x**2+y**2+z**2)
        crot = cos(rot)
        srot = sin(rot)
        costh = abs(crot*z+srot*x)/r
        if(costh.gt.(nt*dcost)) return
        sinth = sqrt(one-costh**2)
        dcosth = dcost
        cosmx = nt*dcos9
        if(jg.le.ng) then
          if(r.le.rmx) then
            ig = jg
            dr = dr0*twig
            lr = nr
            gwt = one
          elseif(jg.lt.ng) then
            ig = jg+1
            dr = two*dr0*twig
            lr = nr
            gwt = eighth
          else
            return
          endif
        else
          if(r.gt.rmx) then
            ig = 2
            dr = four*dr0
            lr = nr
            gwt = (dz9/(four*dz0))*(dx9/(four*dx0))**2
          elseif(abs(costh).gt.cosmx) then
            ig = 1
            dr = two*dr0
            lr = nr
            gwt = (dz9/(two*dz0))*(dx9/(two*dx0))**2
          else
            ig = ng9
            dr = dr9
            lr = nr9
!            gwt = (dz9/(two*dz0))*(dx9/(two*dx0))**2
            gwt = one/nn9
            dcosth = dcos9
          endif
        endif
        if(igr.lt.0) gwt = gwt/four

        ri = r/dr+half
        ir = nint(ri)
        if((ir.lt.1).or.(ir.gt.lr)) return
        dir = ri-ir
        if(dir.lt.zero) then
          dir = -dir
          jr = max(ir-1,1)
        else
          jr = min(ir+1,lr)
        endif
        ti = nt-abs(costh)/dcosth+half
        it = nint(ti)
        if((it.lt.1).or.(it.gt.nt)) return
        dit = ti-it
        if(dit.lt.zero) then
          dit = -dit
          jt = max(it-1,1)
        else
          jt = min(it+1,nt)
        endif
        dirt = one-dir-dit
        if(igr.lt.0) gwt = gwt/four
        if(iy.gt.ly) then
          if(abs(x).lt.rstar) then
            gwt = sqrt(one-(x/rstar)**2)
          else
            gwt = zero
          endif
          if(r.lt.rstar) return
        endif
        if((it.eq.kt).and.(ir.eq.kr).and.(ig.eq.kg)) then
          newcell = .false.
        else
          newcell = .true.
          kt = it
          kr = ir
          kg = ig
        endif

!	let star absorb, but not emit
!	put starlight in at start of starlight rays
        if(starlt) then
          continue
        elseif((jg.eq.ng9).and.(z.gt.zero).and.(z.lt.dz)) then
!	crossing midplane
          taucont = zero
          xi = abs(x)-dx2
          yi = abs(y)-dx2
          ri = sqrt(xi**2+yi**2)
          if(ri.lt.rstar)then
!	inner corner passes through star
            xo = xi+dx
            yo = yi+dx
            ro = sqrt(xo**2+yo**2)
            if(ro.lt.rstar) then
!	entire cell passes through star
              astar = one
            elseif((abs(x).lt.dx).and.(abs(y).lt.dx)) then
!	on one of inner 4 cells
              if(rstar.lt.dx) then
                astar = (pi/4.)*rstar**2/dx**2
              else
                astar = one-(sqrt2-rstar/dx)**2
              endif
            else
!	calculate fraction of cell passing through star
              ri = sqrt(x**2+y**2)-dx2
              ro = ri+dx
              if(ri.gt.rstar) then
                astar = zero
              elseif(ro.lt.rstar) then
                astar = one
              else
                astar = (rstar-ri)/dx
              endif
            endif
            if(astar.gt.zero) then
!	part of cell passed through star
              do iu = 1,nu
                fwi(iu) = fwi(iu)*(one-astar)
              enddo
              if((it.lt.1).or.(it.gt.nt)) return
              do iu = 1,nu
                ff(iu,it,ir,ig) = ff(iu,it,ir,ig)+gwt*fwi(iu)
              enddo
              fn(it,ir,ig) = fn(it,ir,ig)+gwt
              return
            endif
          else
            astar = zero
          endif
        elseif(r.lt.rstar) then
!	still in star off of midplane
          do iu = 1,nu
            fwi(iu) = zero
          enddo
          return
        endif

!	not on star
!	radiative transfer in cell
!	optical depth reduction factor is built into opa and sca
      if(opa(iuv,it,ir,ig)*ds.gt.tiny) then
        do iu = 1,nu
          fwi0 = fwi(iu)
          opac = opa(iu,it,ir,ig)
          scat = sca(iu,it,ir,ig)
          if(iy.le.ly) then
            sfun = swd(iu,it,ir,ig)
          else
            sfun = zero
          endif
!	interpolation across sublimation front causes problems
!          if((ig.gt.jg).and.(ig.le.ng)) then
!            opac = dirt*opac+dir*opa(iu,it,jr,ig)+dit*opa(iu,jt,ir,ig)
!            scat = dirt*scat+dir*sca(iu,it,jr,ig)+dit*sca(iu,jt,ir,ig)
!            if(iy.le.ly) &
!              sfun = dirt*sfun+dir*swd(iu,it,jr,ig)+dit*swd(iu,jt,ir,ig)
!          endif
          if(sfun.lt.zilch) sfun = zero
          opacs = opac+scat
          taud = opacs*ds
          if(taud.lt.small) then
            fwi1 = fwi0+(sfun-fwi0)*taud
            fwi2 = (fwi0+fwi1)/two
          else
!            expd = sexp(-taud)
            arg = hun*taud
            if(arg.lt.big) then
              iarg = int(arg)
              darg = (arg-iarg)/hun
              expd = (one-darg+darg**2/too)*expi(iarg+1)
              fwi1 = fwi0*expd+sfun*(one-expd)
            else
              fwi1 = sfun
            endif
!	Yes.  sfun does belong here.
            fwi2 = (fwi0-fwi1)/taud+sfun
          endif
          if(fwi1.lt.zilch) then
            fwi1 = zero
            fwi2 = fwi0/two
          endif
          fwi(iu) = fwi1
          if(starlt) then
            ff(iu,it,ir,ig) = ff(iu,it,ir,ig)+gwt*fwi2
          else
!	try using incoming radiation field?
!          ff(iu,it,ir,ig) = ff(iu,it,ir,ig)+gwt*fwi0
!	try using radiation field entering r-th cell?
            if(newcell) fwk(iu) = fwi0
            ff(iu,it,ir,ig) = ff(iu,it,ir,ig)+gwt*fwk(iu)
          endif
        enddo  !  iu
      else
        do iu = 1,nu
          if(newcell) fwk(iu) = fwi(iu)
          ff(iu,it,ir,ig) = ff(iu,it,ir,ig)+gwt*fwk(iu)
        enddo
      endif

        fn(it,ir,ig) = fn(it,ir,ig)+gwt

        return
      end



      subroutine gastemp(it,ir,ig,iamax,tsubli,tgi,tdi,hcgas)
!	calculate equilibrium gas and dust temperatures
!	holding pressure constant

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /heat/ gmmdot,acom,pxr,axr,bxr,epspe,agas,sigma, &
                colin,colup,xco,xwat
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg), &
                ff(mu,mt,mr,mg),fn(mt,mr,mg)
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu), &
                phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg), &
                opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg), &
                rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml), &
                wl(ml),dl(ml,mt,mr,mg),sl(ml,mt,mr,mg),opl(ml,mt,mr,mg)
        common /taus/ taucont,tauline,taucent,iuv,iusil,iufl,iwterm
        dimension tda(ma),popa(ma),opaci(mu),scati(mu)
        logical doprint
        logical badnum

        parameter (small = 1.0d-04)
        parameter (zilch = 1.0d-100)
        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (quart = 0.25d+00)
        parameter (fifth = 0.2d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (three = 3.0d+00)
        parameter (four = 4.0d+00)
        parameter (five = 5.0d+00)
        parameter (ten = 1.0d+01)
        parameter (hun = 1.0d+02)
        parameter (pi = 3.14159265d+00)
        parameter (bk = 1.3807d-16)
        parameter (sigsb = 5.670d-05)
        parameter (xmg = 4.0d-24)
        tpi = two*pi
        fpi = four*pi
        doprint = (myrank.eq.0)

!	min and max tg, td
        tgmin = ten
        tgmax = tstar
        tdmin = ten
        tdmax = tsubli
        popmin = 1.0d-50

        if(badnum(tdi).or.(tdi.gt.tdmax)) tdi = tdmax
        tdi = max(tdmin,tdi)
        if(badnum(tgi).or.(tgi.gt.tgmax)) tgi = tgmax
        tgi = max(tgmin,tgi)
        tg0 = tgi
        dens0 = dg(it,ir,ig)
        pres = dens0*tgi
        dens = pres/tgi
        if(ig.eq.ng9) then
          dr = dr9
        else
          dr = dr0*two**ig
        endif
        r = dr*(ir+half)
        if((r.le.rstar).or.(dens.eq.zero)) goto 90
!	initialize opac and scat
        dd0 = dd(it,ir,ig)
        popsum = zero
        do ia = 1,na
!	include settling factor in popa
          popa(ia) = condens(ia,it,ir,ig)*popd(ia)*sexp(-dd0*ad(ia)) &
            *pi*ad(ia)**2
          popsum = popsum+popa(ia)
        enddo
        do iu = 1,nu
          opac = zero
          scat = zero
          if(popsum.gt.popmin) then
            do ia = 1,na
              opac = opac+popa(ia)*phab(iu,ia)
              scat = scat+popa(ia)*phsc(iu,ia)
            enddo
          endif
          opaci(iu) = opac
          scati(iu) = scat
          if(nit.eq.-1) then
            opa(iu,it,ir,ig) = dens*opac
            sca(iu,it,ir,ig) = dens*scat
          endif
        enddo
        if((popsum.lt.zilch).or.(dens.lt.small)) then
          tgi = tstar*(rstar/r)**0.4
          tdi = tsubli
          goto 90
        endif
!	heating and cooling powers/gas particle
        hdust = zero
        hpe = zero
        npe = nu-4
        do iu = nu,npe,-1
          opadu = opaci(iu)*dud(iu)
          hdust = hdust+opadu*fpi*gg(iu,it,ir,ig)
        enddo
!	epspe should depend on dens and maybe ad
        if(hdust.gt.zilch) hpe = epspe*hdust
        hdust = hdust-hpe
        do iu = npe-1,1,-1
          opadu = opaci(iu)*dud(iu)
          hdust = hdust+opadu*fpi*gg(iu,it,ir,ig)
        enddo
!	now include ud < ud(1) with modified dud(1)
!        hdust = hdust+opaci(1)*ud(1)*fpi*gg(1,it,ir,ig)/five
!	old x-ray heating formula assuming absorption indep of freq
!        hxr = (pxr/(fpi*r**2))*sexp(-axr*colin)*axr
!	new formula assumes spectrum power = absorption power + 1
        hxr = (pxr/(fpi*r**2))*((axr/two)/(axr*colin/two+one)**2)
        if(sigma.gt.zero) then
          haccr = gmmdot/(fpi*r**3*sigma)
        else
          haccr = zero
        endif
        hgas = hpe+hxr+haccr
        if(hdust.lt.zilch) hdust = zero
        if(hgas.lt.zilch) hgas = zero
        htot = hgas+hdust
        if(htot.le.zero) then
          tgi = tsubli
          tdi = tsubli
          goto 90
        endif
        dens = pres/tgi
        cgas = cooling(dens,tgi,it,ir,ig)
!	I did have a four in here.  I think it was wrong.
!	I''m not sure about the two.
        acgd = dens*popsum*sqrt(two*bk*tgi/xmg)*acom*2.5*bk
        if(acgd.le.zilch) then
!	no dust and/or no gas, or maybe no coupling
          acgd = zero
          if(hgas*cgas.gt.zero) then
!	must have x-ray heating but no dust
            iit = 0
  112       iit = iit+1
              cgas = cooling(dens,tgi,it,ir,ig)
              tgi = tgi*(hgas/cgas)**(half/max(agas,half))
              dens = pres/tgi
            if((iit.lt.16).and.(tgi.gt.tgmin).and.(tgi.lt.tgmax) &
              .and.(abs(hgas/cgas-one).gt.small)) goto 112
            if(cgas.gt.zero) then
              tgi = min(tgmax,max(tgmin,tgi))
            else
              tgi = tgmax
            endif
          else
!	no gas or no gas heating
            hgas = zero
            tgi = tgmin
          endif
          dens = pres/tgi
          if(hdust.gt.zero) then
!	must have dust but no coupling to gas? doesn''t seem likely
            iit = 0
  114       iit = iit+1
!              cdust = opaci(1)*ud(1)*fpi*bw(ud(1),tdi)/four
              cdust = zero
              do iu = 1,nu
                opadu = opaci(iu)*dud(iu)
                cdust = cdust+opadu*fpi*bw(ud(iu),tdi)
              enddo
              tdi = tdi*(hdust/cdust)**fifth
            if((iit.lt.16).and.(abs(hdust/cdust-one).gt.small)) goto 114
          else
!	no dust or no illumination
            hdust = zero
            tdi = tsubli
          endif
        else
!	calculate tg and td including gas-dust heat flow
          iit = 0
  110     iit = iit+1
!            cdust = opaci(1)*ud(1)*fpi*bw(ud(1),tdi)/four
            cdust = zero
            do iu = 1,nu
              opadu = opaci(iu)*dud(iu)
              cdust = cdust+opadu*fpi*bw(ud(iu),tdi)
            enddo
            cgas = cooling(dens,tgi,it,ir,ig)
            dtd = (hdust+acgd*(tgi-tdi)-cdust)/(five*cdust/tdi+acgd)
            agd = max((agas*cgas/tgi+acgd),half*cgas/tgi)
            dtg = (hgas+acgd*(tdi-tgi)-cgas)/agd
            dtmx = tdi/four
            tdi = tdi+min(dtmx,max(-dtmx,dtd))/two
            tdi = min(tdmax,max(tdmin,tdi))
            dtmx = tgi/four
            tgi = tgi+min(dtmx,max(-dtmx,dtg))/two
            tgi = min(tgmax,max(tgmin,tgi))
            dens = pres/tgi
            if(badnum(tdi).or.badnum(tgi).or.(tgi.le.zero)) then
              print '("tgi",4i4,1p,6e10.2/7e10.2)', &
              it,ir,ig,iit,hdust,hgas,cdust,cgas,acgd,agd, &
              dens,popsum,dtd,dtg,tdi,tgi
              tdi = hun
              tgi = hun
              goto 190
            endif
          if((abs(dtg/tgi).gt.small).and.(iit.lt.16)) goto 110
  190   endif
        hcgas = hgas-cgas
        if(badnum(tdi).or.(tdi.gt.tdmax)) tdi = tdmax
        tdi = max(tdmin,tdi)
        if(badnum(tgi).or.(tgi.gt.tgmax)) tgi = tgmax
        tgi = max(tgmin,tgi)
        dens = pres/tgi
!	tg(it,ir,ig) = tgi  !  in main
!	dg(it,ir,ig) = dens  !  in hydro equil calc
!	calculate individual grain temperatures
!	using heat flows per grain cross section (pi*a**2)
        acgd = dens*sqrt(two*bk*tgi/xmg)*acom*2.5*bk
        iamax = 0
        do ia = 1,na
          if(popa(ia).gt.zero) then
            tdia = tdi
            hdia = zero
            do iu = 1,nu
              hdia = hdia+phab(iu,ia)*fpi*gg(iu,it,ir,ig)*dud(iu)
            enddo
            iit = 0
  130       iit = iit+1
              cdia = zero
              do iu = 1,nu
                cdia = cdia+phab(iu,ia)*fpi*bw(ud(iu),tdia)*dud(iu)
              enddo
              dcdia = five*cdia/tdia+acgd
              dtdia = (hdia-cdia+acgd*(tgi-tdia))/dcdia
              dtmx = tdia/four
              tdia = tdia+min(dtmx,max(-dtmx,dtdia))
              if(badnum(tdia).or.(tdia.le.zero)) then
                print '("tdia",5i4,1p,4e10.2/7e10.2)', &
                  it,ir,ig,ia,iit,hdia,cdia,dcdia,tdia, &
                  tgi,tdi,dens,acgd,hgas,cgas
                stop
              endif
            if((abs(dtdia/tdia).gt.small).and.(iit.lt.16)) goto 130
!	condens is now set in main
!	but should I let small grains sublimate here?
            if(tdia.gt.tsubli) then
              iamax = ia
!              tdia = tsubli
            endif
            tda(ia) = tdia
          else
            tda(ia) = tdi
          endif
        enddo

        if(nit.lt.0) then
          if((taucont.eq.zero).and.(popsum.gt.zero)) then
!	crossing sublimation front; set taufac
            opacv = opaci(iuv)
            scatv = scati(iuv)
            dtauv = sqrt(opacv*(opacv+scatv))*dens*dr
            krmx = 16
            dtaumx = one/krmx
            if(ig.ne.ng9) krmx = krmx/nn9
            if(dtauv.gt.dtaumx) then
              taufac = dtaumx/dtauv
              dtaufac = sqrt(two)
              taufac = taufac/dtaufac
              krmx = min(krmx,int(-log(taufac)/log(dtaufac)))
              do kr = 1,krmx
                taufac = taufac*dtaufac
                jr = ir+kr-1
                if(ig.eq.ng9) then
                  if(jr.le.nr9) then
                    opl(ml,it,jr,ig) = taufac
                  else
                    jr = (jr-1)/(2*nn9)+1
                    jt = nt-(nt-it)/nn9
                    opl(ml,jt,jr,2) = taufac
                  endif
                else
                  if(jr.le.nr) then
                    opl(ml,it,jr,ig) = taufac
                  elseif(ig.lt.ng) then
                    jg = ig+1
                    jr = (jr-1)/2+1
                    opl(ml,it,jr,jg) = taufac
                  endif
                endif
              enddo
            endif
          endif
        endif
        densfac = dens*opl(ml,it,ir,ig)
        dd0 = dd(it,ir,ig)
        do iu = 1,nu
          opac = zero
          scat = aero
          sfun = zero
          do ia = 1,na
            popai = condens(ia,it,ir,ig)*popd(ia)*sexp(-dd0*ad(ia)) &
                    *pi*ad(ia)**2
            if(popai.gt.zero) then
              opac = opac+popai*phab(iu,ia)
              scat = scat+popai*phsc(iu,ia)
              sfun = sfun+popai*phab(iu,ia)*bw(ud(iu),tda(ia))
            endif
          enddo
          opa(iu,it,ir,ig) = opac*densfac
          sca(iu,it,ir,ig) = scat*densfac
          if((opac+scat).gt.zero) then
            swd(iu,it,ir,ig) = &
              (sfun+scat*gg(iu,it,ir,ig))/(opac+scat)
          else
            swd(iu,it,ir,ig) = gg(iu,it,ir,ig)
          endif
        enddo

!	tgi was oscillating.  this seems to cure it.
        tgi = (tg0+tgi)/two
        return

   90   tdi = min(tdmax,max(tdmin,tdi))
        tgi = min(tgmax,max(tgmin,tgi))
        hcgas = zero
        do iu = 1,nu
          w = ud(iu)
          opa(iu,it,ir,ig) = zero
          sca(iu,it,ir,ig) = zero
          swd(iu,it,ir,ig) = bw(w,tdi)
        enddo
        iamax = na+1

        return

      end



      function cooling(dg0,tg0,it,ir,ig)
!	gas cooling by H2O, CO, and H2 Neufeld & Kaufman, 1993
!	per gas particle, T dependance at constant pressure

        implicit real*8 (a-h,o-z)
        parameter (mt = 192)
        parameter (mr = 128)
        parameter (mg = 20)
        parameter (ma = 32)
        parameter (mu = 64)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /parms/ dx0,dz0,dr0,dw,dvturb0,dvturb1,vterm,rstar,tstar
        common /heat/ gmmdot,acom,pxr,axr,bxr,epspe,agas,sigma, &
                colin,colup,xco,xwat
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg), &
                ff(mu,mt,mr,mg),fn(mt,mr,mg)
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu), &
                phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg), &
                opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg), &
                rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (four = 4.0d+00)
        parameter (ten = 1.0d+01)
        parameter (pi = 3.14159265d+00)
        parameter (bk = 1.3807d-16)
        parameter (hck = 1.439d+00)
        parameter (au = 1.496d+13)
        logical badnum
        tpi = two*pi
        fpi = four*pi

        pres = dg0*tg0
        if(badnum(pres).or.(pres.le.zero)) then
          cooling = zero
          agas = one
          return
        endif
        if(colin*colup.gt.zero) then
          colg = one/(one/colin+one/colup)
        else
          colg = zero
        endif
        ds = colg/dg0
        deltat = 1.0d-03*tg0
        call startgr(ig,jg,kr,lr,lt,dr,dcosth)
        r = dr*(ir-half)
        dvturb = sqrt(dvturb0**2+dvturb1**2*au/r)

!	H2O rotational cooling
        if(xwat.gt.zero) then
          tgref = 4.0d+02
          colref = 1.0d+24/dvturb
          tgas = tg0
          dens = pres/tgas
          xmol = xwat
          dmol = xmol*dens
!	c0 = ten**(-23.42)*(tgi/400.)**((23.87-22.88)/log10(1000/200))
          tgif = log10(tgas/tgref)
          cold = log10(xmol*colg/colref+1.0d-04)
          c0 = ten**(-23.42+1.42*tgif)
          clte = ten**(-16.07+3.39*tgif-0.78*cold)
          cn = ten**(4.10+2.19*tgif-cold)
          ca = 0.39-0.20*tgif+0.01*cold
!	assume dust is much less opaque than lines
          cwat = dmol/(one/c0+dens/clte+(dens/cn)**ca*(one/c0-cn/clte))
!	repeat with 1.001*tgas for T dependance
          tgas = tg0+deltat
          dens = pres/tgas
          dmol = xmol*dens
          tgif = log10(tgas/tgref)
!	cold = log10(xmol*colg/colref)
          c0 = ten**(-23.42+1.42*tgif)
          clte = ten**(-16.07+3.39*tgif-0.78*cold)
          cn = ten**(4.10+2.19*tgif-cold)
          ca = 0.39-0.20*tgif+0.01*cold
          dwat = dmol/(one/c0+dens/clte+(dens/cn)**ca*(one/c0-cn/clte))
        else
          cwat = zero
          dwat = zero
        endif
        if(badnum(cwat).or.(cwat.lt.zero)) &
        print '("cwat",3i4,1p,5e10.2/5e10.2)', &
        it,ir,ig,dens,tgas,xmol,tgif,cold,c0,clte,cn,ca,cwat

!	CO rotational cooling
        if(xco.gt.zero) then
          tgref = 6.0d+02
          colref = 1.0d+24/dvturb
          tgas = tg0
          dens = pres/tgas
          xmol = xco
          dmol = xmol*dens
          tgif = log10(tgas/tgref)
          cold = log10(xmol*colg/colref+1.0d-02)
          c0 = ten**(-23.07+0.844*tgif)
          clte = ten**(-18.62+2.72*tgif-0.69*cold)
          cn = ten**(2.79+1.91*tgif-cold)
          ca = 0.46-0.08*tgif+0.06*cold
          cco = dmol/(one/c0+dens/clte+(dens/cn)**ca*(one/c0-cn/clte))
!	repeat with 1.001*tgas for T dependance
          tgas = tg0+deltat
          dens = pres/tgas
          dmol = xmol*dens
          tgif = log10(tgas/tgref)
          c0 = ten**(-23.07+0.844*tgif)
          clte = ten**(-18.62+2.72*tgif-0.69*cold)
          cn = ten**(2.79+1.91*tgif-cold)
          ca = 0.46-0.08*tgif+0.06*cold
          dco = dmol/(one/c0+dens/clte+(dens/cn)**ca*(one/c0-cn/clte))
        else
          cco = zero
          dco = zero
        endif
        if(badnum(cco).or.(cco.lt.zero)) &
        print '("cco",3i4,1p,5e10.2/5e10.2)', &
        it,ir,ig,dens,tgas,xmol,tgif,cold,c0,clte,cn,ca,cco

!	H2 rotational cooling
          tgref = 5.0d+02
          tgas = tg0
          dens = pres/tgas
          dmol = dens
!	just a guess
          umol = 2000.
          iumol = 16.*(slog10(umol)-one)
          taudmol = sexp(-opa(iumol,it,ir,ig)*ds)
          tgif = log10(tgas/tgref)
          c0 = ten**(-24.64+2.27*tgif)
          clte = ten**(-22.13+3.85*tgif)
          cn = ten**(1.14+1.05*tgif)
          ca = 0.42-0.03*tgif
          chh = dmol*taudmol &
            /(one/c0+dens/clte+(dens/cn)**ca*(one/c0-cn/clte))
!	repeat with 1.001*tgas for T dependance
          tgas = tg0+deltat
          dens = pres/tgas
          dmol = dens
          tgif = log10(tgas/tgref)
          c0 = ten**(-24.64+2.27*tgif)
          clte = ten**(-22.13+3.85*tgif)
          cn = ten**(1.14+1.05*tgif)
          ca = 0.42-0.03*tgif
          dhh = dmol*taudmol &
            /(one/c0+dens/clte+(dens/cn)**ca*(one/c0-cn/clte))
        if(badnum(chh).or.(chh.lt.zero)) &
        print '("chh",3i4,1p,6e10.2)', it,ir,ig,dens,c0,clte,cn,ca,chh

!	H2O vibrational cooling
        if(xwat.gt.zero) then
          tgref = 4.0d+02
          colref = 1.0d+24/dvturb
          tgas = tg0
          dens = pres/tgas
          xmol = xwat
          dmol = xmol*dens
          umol = 1617.
          iumol = 16.*(slog10(umol)-one)
          taudmol = sexp(-opa(iumol,it,ir,ig)*ds)
          tgif = log10(tgas/tgref)
          cold = log10(xmol*colg/colref+1.0d-04)
          c0 = 1.03d-26*tgas*sexp(-47.5/tgas**0.333-2325./tgas)
          clte = ten**(-12.99+1.36*tgif-0.74*cold)*exp(-2325./tgas)
          ewat = dmol*taudmol/(one/c0+dens/clte)
!	repeat with 1.001*tgas for T dependance
          tgas = tg0+deltat
          dens = pres/tgas
          dmol = xmol*dens
          tgif = log10(tgas/tgref)
          c0 = 1.03d-26*tgas*sexp(-47.5/tgas**0.333-2325./tgas)
          clte = ten**(-12.99+1.36*tgif-0.74*cold)*exp(-2325./tgas)
          fwat = dmol*taudmol/(one/c0+dens/clte)
        else
          ewat = zero
          fwat = zero
        endif
        if(badnum(fwat).or.(fwat.lt.zero)) &
        print '("fwat",3i4,1p,5e10.2)', it,ir,ig,dens,c0,clte,fwat

!	CO vibrational cooling
        if(xco.gt.zero) then
          tgref = 4.0d+02
          colref = 1.0d+24/dvturb
          tgas = tg0
          dens = pres/tgas
          xmol = xco
          dmol = xmol*dens
          umol = 2142.
          iumol = 16.*(slog10(umol)-one)
          taudmol = sexp(-opa(iumol,it,ir,ig)*ds)
          tgif = log10(tgas/tgref)
          cold = log10(xmol*colg/colref+1.0d-04)
          c0 = 1.83d-26*tgas*sexp(-68.0/tgas**0.333-3080./tgas)
          clte = ten**(-13.44+0.96*tgif-0.83*cold)*exp(-3080./tgas)
          eco = dmol*taudmol/(one/c0+dens/clte)
!	repeat with 1.001*tgas for T dependance
          tgas = tg0+deltat
          dens = pres/tgas
          dmol = xmol*dens
          tgif = log10(tgas/tgref)
          c0 = 1.83d-26*tgas*sexp(-68.0/tgas**0.333-3080./tgas)
          clte = ten**(-13.44+0.96*tgif-0.83*cold)*exp(-3080./tgas)
          fco = dmol*taudmol/(one/c0+dens/clte)
        else
          eco = zero
          fco = zero
        endif
        if(badnum(fco).or.(fco.lt.zero)) &
        print '("fco",3i4,1p,5e10.2)', it,ir,ig,dens,c0,clte,fco

!	H2 vibrational cooling
          tgas = tg0
          dens = pres/tgas
          dmol = dens
          umol = 4680.
          iumol = 16.*(slog10(umol)-one)
          taudmol = sexp(-opa(iumol,it,ir,ig)*ds)
          tgif = log10(tgas/tgref)
          c0 = 1.19d-24*sqrt(tgas) &
            *sexp(-1.81d+04/(tgas+1190.)-5976./tgas)
          clte = 1.10d-18*exp(-6744./tgas)
          ehh = dmol*taudmol/(one/c0+dens/clte)
!	repeat with 1.001*tgas for T dependance
          tgas = tg0+deltat
          dens = pres/tgas
          dmol = dens
          c0 = 1.19d-24*sqrt(tgas) &
            *sexp(-1.81d+04/(tgas+1190.)-5976./tgas)
          clte = 1.10d-18*exp(-6744./tgas)
          fhh = dmol*taudmol/(one/c0+dens/clte)
        if(badnum(fhh).or.(fhh.lt.zero)) &
        print '("fhh",3i4,1p,5e10.2)', it,ir,ig,dens,c0,clte,fhh

        cooling = cwat+cco+chh+ewat+eco+ehh
        dool = dwat+dco+dhh+fwat+fco+fhh
        if(cooling*deltat.gt.zero) then
          agas = (tg0/cooling)*(dool-cooling)/deltat
        else
          agas = one
        endif
        if(badnum(cooling).or.(cooling.lt.zero) &
          .or.badnum(agas).or.(agas.lt.zero)) then
          print '(3i4,1p,2e9.2/8e9.2)', &
            it,ir,ig,tg0,dg0,cooling,agas,cwat,cco,chh,ewat,eco,ehh
          stop
        endif

        return
      end



      subroutine ggavg(nul,ntsm,nrsm)

        implicit real*8 (a-h,o-z)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mg = 20)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /angls/ dcost,dcos9,rot,irot
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg), &
                ff(mu,mt,mr,mg),fn(mt,mr,mg)
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu), &
                phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg), &
                opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg), &
                rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml), &
                wl(ml),dl(ml,mt,mr,mg),sl(ml,mt,mr,mg),opl(ml,mt,mr,mg)
        dimension ggsm(mu,mt),gnsm(mt)
        logical doprint,even
        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (four = 4.0d+00)

        doprint = (myrank.eq.0)
        nr2 = nr/2

        if(ntsm.gt.1) then
          nsm2 = (ntsm-1)/2
          even = (nsm2.ne.(ntsm/2))
!	include angular smoothing in grid 9, but theta_sm will be smaller
          do igr = 0,ng
!	do only 2 point smooth in grids 1-8
            if(igr.gt.0) then
              nsm2 = 0
              even = .true.
            endif
            call startgr(igr,ig,kr,lr,lt,dr,dcosth)
            xlr = lr
            do ir = 1,lr
              do it = 1,lt
                jt1 = it-nsm2
                jt2 = it+nsm2
                gnsm(it) = zero
                do iu = 1,nul
                  ggsm(iu,it) = zero
                enddo
                do jt = jt1,jt2
                  if(jt.lt.1) then
                    kt = 1
                  elseif(jt.gt.nt) then
                    kt = 2*nt-jt+1
                  elseif(jt.gt.lt) then
                    kt = lt
                  else
                    kt = jt
                  endif
                  gnsm(it) = gnsm(it)+gn(kt,ir,ig)
                  do iu = 1,nul
                    ggsm(iu,it) = ggsm(iu,it)+gg(iu,kt,ir,ig)
                  enddo
                enddo
                if(even) then
                  jt1 = jt1-1
                  jt2 = jt2+1
                  if(jt1.lt.1) then
                    jt1 = 1
                  elseif(jt2.gt.nt) then
                    jt2 = 2*nt-jt2+1
                  elseif(jt2.gt.lt) then
                    jt2 = lt
                  endif
                  gnsm(it) = gnsm(it)+(gn(jt1,ir,ig)+gn(jt2,ir,ig))/two
                  do iu = 1,nul
                    ggsm(iu,it) = ggsm(iu,it) &
                      +(gg(iu,jt1,ir,ig)+gg(iu,jt2,ir,ig))/two
                  enddo
                endif
              enddo  !  it
              rrat = (xlr/ir)**2
              do it = 1,lt
                gn(it,ir,ig) = rrat*gnsm(it)
                do iu = 1,nul
                  gg(iu,it,ir,ig) = rrat*ggsm(iu,it)
                enddo
              enddo  !  it
            enddo  !  ir
          enddo  !  igr
        endif

!	smooth mean intensities radially 2 points
        if(nrsm.gt.1) then
          xn9 = two/(one+two*nn9)
          kr2 = nr/2+1
          do igr = 0,ng
            call startgr(igr,ig,kr,lr,lt,dr,dcosth)
            kr = max(kr,2)
            do it = 1,lt
              it1 = nt-(nt-it)/nn9
              do iu = 1,nul
                ggj = gg(iu,it,kr-1,ig)
                ggk = gg(iu,it,kr,ig)
                do ir = kr,lr-1
                  ggi = ggj
                  ggj = ggk
                  ggk = gg(iu,it,ir+1,ig)
                  gg(iu,it,ir,ig) = (ggi+two*ggj+ggk)/four
                enddo
                if(ig.lt.ng) then
                  gg(iu,it,lr,ig) = ((ggi+ggj)/two+two*ggk &
                    +gg(iu,it,kr2,ig+1))/four
                elseif(ig.eq.ng9) then
                  ggi = gg(iu,it1,kr2,2)
                  gg(iu,it,lr,ig) = (ggj+(3.-xn9)*ggk+xn9*ggi)/four
                endif
              enddo
              gnj = gn(it,kr-1,ig)
              gnk = gn(it,kr,ig)
              do ir = kr,lr-1
                gni = gnj
                gnj = gnk
                gnk = gn(it,ir+1,ig)
                gn(it,ir,ig) = (gni+two*gnj+gnk)/four
              enddo
              if(ig.lt.ng) then
                gn(it,lr,ig) = ((gni+gnj)/two+two*gnk &
                  +gn(it,kr2,ig+1))/four
              elseif(ig.eq.ng9) then
                gni = gn(it1,kr2,2)
                gn(it,lr,ig) = (gnj+(3.-xn9)*gnk+xn9*gni)/four
              endif
            enddo  !  it
          enddo  !  igr
        endif

!	divide gg by gn
        gnmin = 1.0d+12
        gnmax = zero
        gnsum = zero
        sumi = zero
        nzero = 0
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          ls = nt-lt+1
          do ir = kr,lr
            r = dr*(ir-half)
            do it = 1,lt
              gni = gn(it,ir,ig)
              if(gni.eq.zero) then
                if((r.gt.rstar).and.(opl(1,it,ir,ig).gt.zero)) then
                  if(doprint.and.(nzero.eq.0)) &
                    print'("Warning: no rays pass through ig,ir,it =", &
                    3i6)', ig,ir,it
                  nzero = nzero+1
                endif
              elseif((ir.gt.lr/2).and.(it.gt.(nt/2))) then
                if(gni.lt.gnmin) then
                  gnmin = gni
                  itmin = it
                  irmin = ir
                  igmin = ig
                elseif(gni.gt.gnmax) then
                  gnmax = gni
                  itmax = it
                  irmax = ir
                  igmax = ig
                endif
                gnsum = gnsum+gni
                sumi = sumi+one
              endif
              if(gni.gt.zero) then
                do iu = 1,nul
                  gg(iu,it,ir,ig) = gg(iu,it,ir,ig)/gni
                enddo
              elseif(ir.gt.kr) then
                do iu = 1,nul
                  gg(iu,it,ir,ig) = gg(iu,it,ir-1,ig)
                enddo
              else
                do iu = 1,nul
                  gg(iu,it,ir,ig) = zero
                enddo
              endif
            enddo  !  it
          enddo  !  ir
        enddo  !  igr
        gnmean = gnsum/sumi
        if(gnmin.gt.gnmax) gnmin = zero
        if(doprint.and.(nal.eq.na)) then
          print '(" min,max,mean rays per cell:",1p,3f12.2)', &
            gnmin,gnmax,gnmean
          print '(6i8)', itmin,irmin,igmin,itmax,irmax,igmax
          if(nzero.gt.0) &
          print '(" no rays pass through",i8," cells of",i8)', &
            nzero,nt*nr2*(ng-1)+nt*nr9
        endif

        return
      end



      subroutine initju(ju,kk)
!	initialize variables dependent on ju

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        parameter (mw = 1024)
        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /xj/ xju,xjp,xjq,xjr,xk
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg), &
                ff(mu,mt,mr,mg),fn(mt,mr,mg)
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu), &
                phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg), &
                opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg), &
                rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml), &
                wl(ml),dl(ml,mt,mr,mg),sl(ml,mt,mr,mg),opl(ml,mt,mr,mg)
        logical doprint
        parameter (half = 0.5d+00)
        parameter (two = 2.0d+00)

        doprint = (myrank.eq.0)

        if(moltype.eq.3) then
!	for NH3
          if(mod(kk,3).eq.0) then
            gnucl = 2.
          else
            gnucl = 1.
          endif
        elseif(moltype.eq.4) then
!	for C2H2 v4+v5
          if(mod(ju,2).eq.0) then
            gnucl = 0.5
          else
            gnucl = 1.5
          endif
        elseif(moltype.eq.5) then
!	for C2H2 P+R
!	opacities can be corrected to 2-gnucl for Q in step.  or not
          if(mod(ju,2).eq.0) then
            gnucl = 1.5
          else
            gnucl = 0.5
          endif
        else
          gnucl = 1.
        endif
        jr = ju-1
        xju = ju
        xjr = jr
        wr0 = two*xju*be
        wl(1) = wr0

!	initialize dust sfun, opac, molec densities
!	interpolate from swd and opa
!	should I do all lines here?
        wi = wl(1)
!	w = w0*10.**(dlogw*iu)
        ui = slog(wi/ud(1))/slog(ud(2)/ud(1))+1.
        iu1 = ui-1.
        iu2 = iu1+1
        iu3 = iu1+2
        iu4 = iu1+3
        if(iu1.lt.1) then
          if(doprint) print &
            '("Line outside of continuum range.  Try extrapolating")'
          iu1 = 1
          iu2 = 2
          iu3 = 3
          iu4 = 4
          du1 = ui-iu1
          du2 = ui-iu2
          wt1 = -du2
          wt2 = du1
          wt3 = zero
          wt4 = zero
        elseif(iu4.gt.nu) then
          if(doprint) then
            print '("Line outside of continuum range")'
            stop
          else
            call sleep(10)
          endif
        else
          du1 = ui-iu1
          du2 = ui-iu2
          du3 = ui-iu3
          du4 = ui-iu4
          wt1 = -du2*du3*du4/6.
          wt2 = du1*du3*du4/2.
          wt3 = -du1*du2*du4/2.
          wt4 = du1*du2*du3/6.
        endif
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            do it = 1,lt
              sw1 = slog(swd(iu1,it,ir,ig))
              sw2 = slog(swd(iu2,it,ir,ig))
              sw3 = slog(swd(iu3,it,ir,ig))
              sw4 = slog(swd(iu4,it,ir,ig))
              swi = sw1*wt1+sw2*wt2+sw3*wt3+sw4*wt4
              sl(1,it,ir,ig) = sexp(swi)
!	should interpolate gg for rotational lines if I want to include
!	radiative excitation
              op1 = opa(iu1,it,ir,ig)
              op2 = opa(iu2,it,ir,ig)
              op3 = opa(iu3,it,ir,ig)
              op4 = opa(iu4,it,ir,ig)
              opl(1,it,ir,ig) = op1*wt1+op2*wt2+op3*wt3+op4*wt4
!	neglect scattering by dust at the line wavelengths
            enddo
          enddo
        enddo

        return
      end



      subroutine grid0(fwk,ly)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwk(mw)
        dimension fwi(mw)
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /angls/ dcost,dcos9,rot,irot
        common /taus/ taucont,tauline,taucent,iuv,iusil,iufl,iwterm
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

!	first grid to calculate
        ig = ng
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

!	zero spectrum summed over lines of sight
        do iw = 1,lw
          fwk(iw) = zero
        enddo
        do ix = 1,nx
          do iy = 1,ly
!	in middle 1/2 x 1/2 of grid?
            middle = (ix.gt.nx4).and.(ix.le.nx3) &
              .and.(iy.gt.nx4).and.(iy.le.nx3)
            if(middle) then
!	use finer grid for grid section above and below middle
              jx2 = 2*ix-nx2
              jx1 = jx2-1
              jy2 = 2*iy-nx2
              jy1 = jy2-1
              do jx = jx1,jx2
                do jy = jy1,jy2
!	zero spectrum for line of sight
                  do iw = 1,lw
                    fwi(iw) = zero
                  enddo
                  iz = 0
                  wdopp = zero
!	step along line of sight
  14              iz = iz+1
                  if(iz.eq.nz5) then
                    if(ig.eq.12) then
                      call gridb(fwi,jx,jy,wdopp)
                    elseif(ig.eq.11) then
                      call grida(fwi,jx,jy,wdopp)
                    elseif(ig.eq.10) then
                      call grid9(fwi,jx,jy,wdopp)
                    elseif(ig.eq.9) then
                      call grid8(fwi,jx,jy,wdopp)
                    elseif(ig.eq.8) then
                      call grid7(fwi,jx,jy,wdopp)
                    elseif(ig.eq.7) then
                      call grid6(fwi,jx,jy,wdopp)
                    elseif(ig.eq.6) then
                      call grid5(fwi,jx,jy,wdopp)
                    elseif(ig.eq.5) then
                      call grid4(fwi,jx,jy,wdopp)
                    elseif(ig.eq.4) then
                      call grid3(fwi,jx,jy,wdopp)
                    elseif(ig.eq.3) then
                      call grid2(fwi,jx,jy,wdopp)
                    else
                      print '(" ng =",i4,"?")', ig
                      stop
                    endif
                    iz = iz+nz/2
                  endif
!	call radiative transfer routine for step
                  call step(fwi,wdopp,jx,jy,iz,-ig)
!	end loop along line of sight
                  if(iz.lt.nz) go to 14
!	store compressed spectrum in data cube
                  if(ig.le.ngc) call storecube(fwi,jg,jx,jy)
!	add outgoing spectrum to sum from grid
                  do iw = 1,lw
                    fwk(iw) = fwk(iw)+fwi(iw)/4.
                  enddo
!	end loops over 4 lines of sight
                enddo
              enddo
            else
!	not in middle.  follow ray along ix,iy
              do iw = 1,lw
                fwi(iw) = zero
              enddo
              wdopp = zero
!	step along line of sight
              do iz = 1,nz
!	call radiative transfer routine for step
                call step(fwi,wdopp,ix,iy,iz,ig)
              enddo
!	store compressed spectrum in data cube
              if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
!	add outgoing spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwi(iw)
              enddo
            endif
!	end loops over lines of sight in ix,iy grid
          enddo
        enddo

        return
      end


      subroutine gridb(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 11
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3) &
          .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
!	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
!	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
!	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grida(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
!	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
!	end loop along line of sight
              if(iz.lt.nz) go to 24
!	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
!	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
!	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
!	don''t branch out
          wdopp = wdopp0
!	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
!	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grida(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 10
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3) &
          .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
!	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
!	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
!	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid9(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
!	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
!	end loop along line of sight
              if(iz.lt.nz) go to 24
!	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
!	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
!	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
!	don''t branch out
          wdopp = wdopp0
!	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
!	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid9(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 9
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3) &
          .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
!	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
!	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
!	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid8(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
!	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
!	end loop along line of sight
              if(iz.lt.nz) go to 24
!	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
!	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
!	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
!	don''t branch out
          wdopp = wdopp0
!	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
!	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid8(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 8
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3) &
          .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
!	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
!	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
!	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid7(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
!	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
!	end loop along line of sight
              if(iz.lt.nz) go to 24
!	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
!	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
!	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
!	don''t branch out
          wdopp = wdopp0
!	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
!	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid7(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 7
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3) &
          .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
!	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
!	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
!	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid6(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
!	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
!	end loop along line of sight
              if(iz.lt.nz) go to 24
!	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
!	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
!	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
!	don''t branch out
          wdopp = wdopp0
!	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
!	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid6(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 6
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3) &
          .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
!	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
!	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
!	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid5(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
!	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
!	end loop along line of sight
              if(iz.lt.nz) go to 24
!	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
!	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
!	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
!	don''t branch out
          wdopp = wdopp0
!	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
!	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid5(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 5
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3) &
          .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
!	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
!	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
!	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid4(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
!	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
!	end loop along line of sight
              if(iz.lt.nz) go to 24
!	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
!	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
!	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
!	don''t branch out
          wdopp = wdopp0
!	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
!	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid4(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 4
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3) &
          .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
!	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
!	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
!	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid3(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
!	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
!	end loop along line of sight
              if(iz.lt.nz) go to 24
!	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
!	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
!	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
!	don''t branch out
          wdopp = wdopp0
!	step along line of sight
          do iz = 1,nz
!	call radiative transfer routine for step
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
!	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid3(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 3
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3) &
          .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
!	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
!	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
!	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid2(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
!	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
!	end loop along line of sight
              if(iz.lt.nz) go to 24
!	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
!	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
!	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
!	don''t branch out
          wdopp = wdopp0
!	step along line of sight
          do iz = 1,nz
!	call radiative transfer routine for step
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
!	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid2(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 2
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3) &
          .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
!	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
!	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
!	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid1(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
!	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
!	end loop along line of sight
              if(iz.lt.nz) go to 24
!	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
!	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
!	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
!	not in the middle, don''t branch out
          wdopp = wdopp0
!	step along line of sight
          do iz = 1,nz
!	call radiative transfer routine for step
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
!	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid1(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mu = 64)
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /angls/ dcost,dcos9,rot,irot
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = ng9
!	middle is not a special case for grid 1
!	store expanded dust grid to cube(,,,2)

        do iw = 1,lw
          fwk(iw) = zero
        enddo
!	branch out by nn9*nn9
        jx1 = nn9*(ix-1)+1
        jx2 = nn9*ix
        jy1 = nn9*(iy-1)+1
        jy2 = nn9*iy
        do jx = jx1,jx2
          do jy = jy1,jy2
!	copy spectrum from outer grid
            do iw = 1,lw
              fwj(iw) = fwi(iw)
            enddo
            wdopp = wdopp0
!	step along line of sight
            do iz = 1,nz9
              call step(fwj,wdopp,jx,jy,iz,ig)
            enddo
!	store compressed dust spectrum in data cube
            if(moltype.le.0) call storecube(fwj,ig,jx,jy)
!	add spectrum to sum from grid
            do iw = 1,lw
              fwk(iw) = fwk(iw)+fwj(iw)
            enddo
!	end loops over lines of sight on grid
          enddo
        enddo
!	put sum from grid into fwi
        xn9 = nn9**2
        do iw = 1,lw
          fwi(iw) = fwk(iw)/xn9
        enddo
        wdopp0 = wdopp

        return
        end



      subroutine storecube(fwi,ig,ix,iy)
!	store compressed spectrum in data cube

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        dimension fwi(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /taus/ taucont,tauline,taucent,iuv,iusil,iufl,iwterm
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        if(ig.eq.ng9) then
          if(moltype.gt.0) return
          nc = mc
          if(nc.ge.nx9) then
            ncnx = nc/nx9
            kx4 = ncnx*ix
            kx3 = kx4-ncnx+1
            ky4 = ncnx*iy
            ky3 = ky4-ncnx+1
            if((kx3.lt.1).or.(kx4.gt.mc).or.(ky3.lt.1).or.(ky4.gt.mc)) then
              print '(" in storecube kx3,kx4,ky3,ky4 =",5i8)', &
                kx3,kx4,ky3,ky4,mc
              return
            endif
            fwisil = fwi(iusil)
            do ky = ky3,ky4
              ky1 = nc+1-ky
              do kx = kx3,kx4
                cube(kx,ky,2) = fwisil
                cube(kx,ky1,2) = fwisil
              enddo
            enddo
          else
            ncnx = nx9/nc
            kx = (ix+1)/ncnx
            ky = (iy+1)/ncnx
            if((kx.lt.1).or.(kx.gt.nc).or.(ky.lt.1).or.(ky.gt.nc)) then
              print '(" in storecube kx3,kx4,ky3,ky4 =",5i8)', &
                kx3,kx4,ky3,ky4,nc
              return
            endif
            fwisil = fwi(iusil)/ncnx
            ky1 = nc+1-ky
            if((mod(ix,ncnx).eq.1).and.(mod(iy,ncnx).eq.1)) then
              cube(kx,ky,2) = fwisil
              cube(kx,ky1,2) = fwisil
            else
              cube(kx,ky,2) = cube(kx,ky,2)+fwisil
              cube(kx,ky1,2) = cube(kx,ky1,2)+fwisil
            endif
          endif
        elseif(ig.gt.ngc) then
          return
        else
          nc = mc/2**(ngc-ig)
          ncnx = nc/nx
          middle = (ix.gt.nx/4).and.(ix.le.3*nx/4) &
            .and.(iy.gt.nx/4).and.(iy.le.3*nx/4)
          if((middle.and.(ncnx.gt.1).and.(ig.gt.1)) &
            .or.(ncnx.lt.1)) return
          kx4 = ncnx*ix+(mc-nc)/2
          kx3 = kx4-ncnx+1
          ky4 = ncnx*iy+(mc-nc)/2
          ky3 = ky4-ncnx+1
          if((kx3.lt.1).or.(kx4.gt.mc).or.(ky3.lt.1).or.(ky4.gt.mc)) then
            print '(" in storecube kx3,kx4,ky3,ky4 =",5i8)', &
              kx3,kx4,ky3,ky4,mc
            stop
          endif
          if(moltype.gt.0) then
            nnw = 2*nw/mc
            xnw = nnw
            iw = 1
            do ky = ky3,ky4
              ky1 = mc+1-ky
              do kx = kx3,kx4
                cube(kx,ky,iw) = fwi(1)
                cube(kx,ky1,iw) = fwi(1)
              enddo
            enddo
            do iw0 = nnw,nw-nnw,nnw
              iw = iw0/nnw+1
              fwii = zero
              fwij = zero
              do jw = iw0+1,iw0+nnw
                jw1 = nw+1-jw
                fwii = fwii+fwi(jw)
                fwij = fwij+fwi(jw1)
              enddo
              do ky = ky3,ky4
                ky1 = mc+1-ky
                do kx = kx3,kx4
                  cube(kx,ky,iw) = fwii/xnw
                  cube(kx,ky1,iw) = fwij/xnw
                enddo
              enddo
            enddo
          else
            fwisil = fwi(iusil)
            do ky = ky3,ky4
              ky1 = mc+1-ky
              do kx = kx3,kx4
                cube(kx,ky,1) = fwisil
                cube(kx,ky1,1) = fwisil
              enddo
            enddo
          endif
        endif

        return
      end



      subroutine step(fwi,wdopp0,ix,iy,iz,igr)

        implicit real*8 (a-h,o-z)
        parameter (ml = 2)
        parameter (mw = 1024)

        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /angls/ dcost,dcos9,rot,irot
        dimension fwi(mw)

        if(moltype.le.0) then
          call stepd(fwi,ix,iy,iz,igr)
        else
          call stepl(fwi,wdopp0,ix,iy,iz,igr)
        endif

        return
      end




      subroutine stepl(fwi,wdopp0,ix,iy,iz,igr)

        implicit real*8 (a-h,o-z)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mg = 20)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        parameter (mw = 1024)

        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg), &
                ff(mu,mt,mr,mg),fn(mt,mr,mg)
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /angls/ dcost,dcos9,rot,irot
        common /xj/ xju,xjp,xjq,xjr,xk
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu), &
                phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg), &
                opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg), &
                rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml), &
                wl(ml),dl(ml,mt,mr,mg),sl(ml,mt,mr,mg),opl(ml,mt,mr,mg)
        common /vels/ wr(mt,mr,mg),wth(mt,mr,mg),wph(mt,mr,mg)
        common /taus/ taucont,tauline,taucent,iuv,iusil,iufl,iwterm
        common /exptab/ expi(10239),too
        dimension fwi(mw)
        logical good,doprint,badnum

        parameter (pi = 3.14159265d+00)
        parameter (sqrtpi = 1.772454d+00)
        parameter (sqrt2 = 1.414214d+00)
        parameter (alog2 = 0.6931472d+00)
        parameter (small = 1.0d-03)
        parameter (tiny = 1.0d-08)
        parameter (zilch = 1.0d-200)
        parameter (hun = 1.0d+02)
        parameter (big = 5.0d+03)
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (four = 4.0d+00)
        parameter (half = 0.5d+00)
        parameter (eighth = 0.125d+00)
        parameter (hc2 = 5.955d-06)
        parameter (au = 1.496d+13)
        parameter (clt = 2.9979d+10)

        doprint = (myrank.eq.0)
        nw2 = nw/2
        grat = (two-gnucl)/gnucl
!	set grat = 1 for a pseudo-line with same gnucl as P,R

        jg = abs(igr)
        if((jg.le.1).or.(jg.eq.ng9)) then
          jg = ng9
          lx = nx9
          dx = dx9
          lz = nz9
          dz = dz9
          rmx = nr9*dr9
        else
          twig = two**jg
          lx = nx
          dx = dx0*twig
          lz = nz
          dz = dz0*twig
          rmx = nr*dr0*twig
        endif
        if(igr.lt.0) dx = dx/two
        lx2 = lx/2
        if((igr.ge.0).and.(ix.eq.lx2).and.(iy.eq.1).and.(iz.eq.1)) then
          taucont = zero
          tauline = zero
          taucent = zero
        endif
        dx2 = dx/two
        dz2 = dz/two
        ds = dz
        x = dx2*(2*ix-lx-1)
        y = dx2*(2*iy-lx-1)
        z = dz2*(2*iz-lz-1)
        r = sqrt(x**2+y**2+z**2)
        crot = cos(rot)
        srot = sin(rot)
        costh = (crot*z+srot*x)/r
        if(abs(costh).gt.(nt*dcost)) return
        sinth = sqrt(one-costh**2)
        dcosth = dcost
        cosmx = nt*dcos9
        ph = atan2(y,(crot*x-srot*z))
        if(jg.le.ng) then
          if(r.le.rmx) then
            ig = jg
            dr = dr0*twig
            lr = nr
            gwt = one
          elseif(jg.lt.ng) then
            ig = jg+1
            dr = two*dr0*twig
            lr = nr
            gwt = eighth
          else
            return
          endif
        else
          if(r.gt.rmx) then
            ig = 2
            dr = four*dr0
            lr = nr
            gwt = (dz9/(four*dz0))*(dx9/(four*dx0))**2
          elseif(abs(costh).gt.cosmx) then
            ig = 1
            dr = two*dr0
            lr = nr
            gwt = (dz9/(two*dz0))*(dx9/(two*dx0))**2
            if(igr.eq.1) jg = 1
          else
            ig = ng9
            dr = dr9
            lr = nr9
!            gwt = (dz9/(two*dz0))*(dx9/(two*dx0))**2
            gwt = one/nn9
            dcosth = dcos9
          endif
        endif
        ir = r/dr+one
        if(igr.lt.0) gwt = gwt/four
        it = nt-nint(abs(costh)/dcosth-half)
        if((it.lt.1).or.(it.gt.nt)) return

!	does ray hit star?
        if((jg.eq.ng9).and.(z.gt.zero).and.(z.lt.dz)) then
          xi = abs(x)-dx2
          yi = abs(y)-dx2
          ri = sqrt(xi**2+yi**2)
          if(ri.lt.rstar) then
            xo = xi+dx
            yo = yi+dx
            ro = sqrt(xo**2+yo**2)
            if(ro.lt.rstar) then
              astar = one
            elseif((abs(x).lt.dx).and.(abs(y).lt.dx)) then
              if(rstar.lt.dx) then
                astar = (pi/4.)*rstar**2/dx**2
              else
                astar = one-(sqrt2-rstar/dx)**2
              endif
            else
              ri = sqrt(x**2+y**2)-dx2
              ro = ri+dx
              if(ri.gt.rstar) then
                astar = zero
              elseif(ro.lt.rstar) then
                astar = one
              else
                astar = (rstar-ri)/dx
              endif
            endif
            if(astar.gt.zero) then
              fstara = fstar(1)*astar
              do iw = 1,nw
                fwi(iw) = (one-astar)*fwi(iw)+fstara
              enddo
!              ff(1,it,ir,ig) = ff(1,it,ir,ig)+gwt*fwi(1)
!              fn(it,ir,ig) = fn(it,ir,ig)+gwt
              return
            endif
          endif
        elseif(r.lt.rstar) then
          do iw = 1,nw
            fwi(iw) = fstar(1)
          enddo
          return
        endif

!	not on star
        wr0 = two*xju*be
        wrc = wr0/clt
        dw = dv*wrc
        dw2 = dw/two
        srd = sl(1,it,ir,ig)
        vr = wr(it,ir,ig)
        vth = wth(it,ir,ig)
        vph = wph(it,ir,ig)
        vtot = sqrt(vr**2+vth**2+vph**2)
        vz = vr*z/r+vph*sin(ph)*srot
        if(vth.ne.zero) vz = vz-vth*(sinth*crot+costh*cos(ph)*srot)
        da = dg(it,ir,ig)
        if(da.lt.zilch) then
          da = zero
          djr = zero
          dju = zero
          rdcoll = zero
        else
          djr = dl(1,it,ir,ig)
          dju = dl(2,it,ir,ig)
        endif

        wdopp = vz*wrc
        dwdopp = abs(wdopp0-wdopp)/two
        wdopp0 = wdopp
        if(iz.eq.1) dwdopp = dwdopp*2./3.
!	I ought to multiply by 4/3 going back out a grid
        dwr = sqrt(dvturb0**2+dvturb1**2*au/r+dv**2)*wrc
        if(dwr.eq.zero) then
          print '(" dvturb = dw = 0.")'
          stop
        endif
        dwhwzi = dwdopp+4.*dwr
        phiint = two*dwdopp+sqrtpi*dwr
        fri = zero
        wr3 = wr0**3
        wmid = (nw+one)/two
        opacdr = opl(1,it,ir,ig)
        taudr = opacdr*ds
        expdr = sexp(-taudr)
        gju = two*xju+one
        gjr = two*xjr+one
        if(moltype.eq.1) then
!	diatomic and linear molecule opacities
          opfacr = (djr/gjr-dju/gju)*xju*sv
          emfacr = two*hc2*wr3*sv*dju*xju/gju
        elseif(moltype.eq.2) then
!	HCN opacities are the same?
          opfacr = (djr/gjr-dju/gju)*xju*sv
          emfacr = two*hc2*wr3*sv*dju*xju/gju
        elseif(moltype.eq.3) then
!	NH3 opacities?
          opfacr = (djr/gjr-dju/gju)*(xju-xk*xk/xju)*sv
          emfacr = two*hc2*wr3*sv*dju*(xju-xk*xk/xju)/gju
        elseif(doprint) then
          print '("Undefined moltype")'
          stop
        endif
        if(emfacr.lt.zilch) emfacr = zero
        good = .true.

        taur0 = opfacr*ds/phiint
        if(taur0.lt.(-tiny)) then
          if(doprint) print '("Warning: tau0 =",1p,e10.2," at",3i6)', &
            taur0,ig,ir,it
        endif
        if(taur0.gt.tiny) then
          dwmax = dwdopp+dwr*sqrt(slog(taur0/tiny))
        else
          dwmax = zero
        endif
        if(dwmax.lt.dw) dwmax = zero
        phisum = zero
!	loop over frequencies
        do iw = 1,nw
          fwr0 = fwi(iw)
          w = (iw-wmid)*dw
          deltaw = abs(w-wdopp)
          if((dwmax.gt.zero).and.(iw.gt.1).and.(iw.lt.nw)) then
!	include line opacity
            if(deltaw.lt.dwmax) then
!	crude convolution of Gaussian with box
              if(deltaw.lt.dwdopp) then
                phiw = one
              else
!	         phiw = sexp(-((deltaw-dwdopp)/dwr)**2)
                arg = hun*((deltaw-dwdopp)/dwr)**2
                if(arg.lt.big) then
                  iarg = int(arg)
                  darg = (arg-iarg)/hun
                  phiw = (one-darg+darg**2/too)*expi(iarg+1)
                else
                  phiw = zero
                endif
              endif
              phiw = phiw/phiint
              phiwdw = phiw*dw
              phisum = phisum+phiwdw
              opacr = opacdr+opfacr*phiw
              emisr = opacdr*srd+emfacr*phiw
              if((emisr.gt.zilch).and.(opacr.gt.zilch)) then
                sr = emisr/opacr
              else
                sr = zero
              endif
              taur = opacr*ds
!	leave out pisco change allowing negative tauline
              if(taur.lt.small) then
                fwr1 = fwr0+(sr-fwr0)*taur
              elseif((taur-taudr).lt.small) then
                expr = expdr*(one-(taur-taudr))
                fwr1 = fwr0*expr+sr*(one-expr)
              else
                arg = hun*taur
                if(arg.lt.big) then
                  iarg = int(arg)
                  darg = (arg-iarg)/hun
                  expr = (one-darg+darg**2/too)*expi(iarg+1)
                  fwr1 = fwr0*expr+sr*(one-expr)
                else
                  expr = zero
                  fwr1 = sr
                endif
              endif
            else
!	include only dust off line
              if(taudr.gt.small) then
                fwr1 = fwr0*expdr+srd*(one-expdr)
              else
                fwr1 = fwr0+(srd-fwr0)*taudr
              endif
              taur = taudr
            endif
          else
!	no line opacity, include only dust
            if(taudr.gt.small) then
              fwr1 = fwr0*expdr+srd*(one-expdr)
            else
              fwr1 = fwr0+(srd-fwr0)*taudr
            endif
            taur = taudr
          endif
!	intensity leaving cell
          fwi(iw) = fwr1
          if(good.and.(badnum(fwr1).or.(fwr1.lt.0.))) then
            print '(" fwr1 bad",2i6,1p,5e10.3)', iw,ig,x,y,z,opacr,emisr
            print '(1p,8e9.2)', taur,taudr,fwr1,fwr0,expr,sr,expdr,srd
            good = .false.
            stop
          endif
!	calculate and print optical depths
          if((igr.ge.0).and.(ix.eq.lx2).and.(iy.eq.1) &
            .and.(deltaw.lt.dw2)) then
            tauline = tauline+taur-taudr
          endif
        enddo  !  iw
        if((igr.ge.0).and.(ix.eq.lx2).and.(iy.eq.1)) then
          taucont = taucont+taudr
          taucent = taucent+opfacr*ds/phiint
          if(doprint.and.(iz.eq.lz)) then
            print '(" at y, rot =",2f7.2, &
              " continuum, line taus =",1p,3e10.2)', &
              -y/au,rot,taucont,tauline,taucent
            write(16,'(" at y, rot =",2f7.2, &
              " continuum, line taus =",1p,3e10.2)') &
              -y/au,rot,taucont,tauline,taucent
          endif
        endif

!	need to calc ff if doing NLTE calculation

        return
      end



      subroutine movie(nc,ic)

        implicit real*8 (a-h,o-z)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /angls/ dcost,dcos9,rot,irot
        common /cube/ cube(mc,mc,mc),ngc
        dimension frame(mc,mc)
        character*16 lockline

!	make and plot continuum frame
        cosrot = cos(rot)
        xmid = (nc+1)/two
        nc2 = nc/2
        sumcon = zero
        do iy = 1,nc
          do ix = 1,nc
            x = (ix-xmid)/cosrot
            y = iy-xmid
            scale = sqrt(x**2+y**2)
            scale = (x**2+y**2)
            if(scale.lt.one) scale = zero
            frame(ix,iy) = sqrt(max(zero,scale*cube(ix,iy,ic)))
            sumcon = sumcon+frame(ix,iy)
          enddo
        enddo
        if(sumcon.eq.zero) then
          print '(" continuum = 0")'
          return
        endif
   18   open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=19)
        read(8,'(a16)') lockline
        close(unit=8)
        if(index(lockline,'reading').gt.0) then
          call sleep(1)
          goto 18
        endif
        open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=19)
        write(8,'("writing")')
        close(unit=8)
        open(unit=9,file='/home/r900-1/lacy/plots/images', &
          access='direct',err=19,recl=2*nc)
        irec = 0
        do iy = 1,nc
          irec = irec+1
          write(unit=9,rec=irec) (frame(ix,iy),ix=1,nc)
        enddo
        nfr = 1
        if(lw.lt.(nw/2)) goto 17

!	make and plot line movie
        nw1 = nc/8
        nw2 = nc/4
        nw3 = 5*nc/16
        idw = nw2/16
        nw1 = nw1-idw/2
        xdw = idw
        do iw = nw2,nw3,idw
          sumsum = zero
          do iy = 1,nc
            do ix = 1,nc
              sumj = -xdw*cube(ix,iy,ic)
              do jw = iw,iw+idw-1
                sumj = sumj+cube(ix,iy,jw)
              enddo
              x = (ix-xmid)/cosrot
              y = iy-xmid
              scale = sqrt(x**2+y**2)/xdw
              scale = (x**2+y**2)/xdw
!              if(scale.lt.one) scale = zero
              frame(ix,iy) = scale*sumj
              sumsum = sumsum+abs(scale*sumj)
            enddo
          enddo
          if(sumsum.gt.(sumcon/1.0d+06)) then
!            print '(" iw =",i4," -",i4,1p,2e10.2)', &
!              iw,iw+idw-1,sumsum,sumcon
            do iy = 1,nc
              irec = irec+1
              write(unit=9,rec=irec) (frame(ix,iy),ix=1,nc)
            enddo
            nfr = nfr+1
          endif
        enddo

!	make line frame
        nc2 = nc/2
        sumsum = zero
        do iy = 1,nc
          do ix = 1,nc
            x = (ix-xmid)/cosrot
            y = iy-xmid
            scale = sqrt(x**2+y**2)
            scale = (x**2+y**2)
!            if(scale.lt.one) scale = zero
            sum = -(nc2-2)*cube(ix,iy,ic)
            do iw = 2,nc2-1
              sum = sum+cube(ix,iy,iw)
            enddo
            sum = scale*sum
            frame(ix,iy) = sum
            sumsum = sumsum+abs(sum)
          enddo
        enddo
        if(sumsum.ne.zero) then
          do iy = 1,nc
            irec = irec+1
            write(unit=9,rec=irec) (frame(ix,iy),ix=1,nc)
          enddo
          nfr = nfr+1
        endif
   17   close(unit=9)
        open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=19)
        write(8,'("images")')
        write(8,'(3i8)') nc,nc,nfr
        close(unit=8)

   19   return
      end



      subroutine convo(fw,spec,fwnorm,dvdpix,wline,dsm,nsm)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        common /mpis/mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        dimension fw(mw),spec(mw)
        logical doprint
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)

        doprint = (myrank.eq.0)
        sumss = zero
        do iw = 1,nw
          sums = fw(iw)
          sumg = one
          do ism = 1,nsm
            iwm = iw-ism
            iwp = iw+ism
            if((iwm.ge.1).and.(iwp.le.nw)) then
              gaus = sexp(-(ism/dsm)**2)
              sumg = sumg+two*gaus
              sums = sums+gaus*(fw(iwm)+fw(iwp))
            endif
          enddo
          speci = sums/(sumg*fwnorm)
          spec(iw) = speci
          sumss = sumss+speci
        enddo
!	fwnorm is the normalization factor to convert fw to Jy
!	or normalize to the stellar continuum
!	if spec is in Jy, fline is in Jy-km/s
        cont = (fw(1)+fw(nw))/(two*fwnorm)
        ew = (nw-sumss/cont)*dvdpix
        fline = -ew*cont
        gline = 1.0d-21*wline*fline
!        print '("      Fc(Jy) Fl(Jy-km/s)    Fl(W/m2)    EW(km/s)")'
        if(doprint) then
          print '(1p,4e12.3)', cont,fline,gline,ew
          write(16,'(1p,4e12.3)') cont,fline,gline,ew
        endif

        return
      end



      subroutine zerof(nul)
!	zero ff and fn

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (ml = 2)
        parameter (mu = 64)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg), &
                ff(mu,mt,mr,mg),fn(mt,mr,mg)
        parameter (zero = 0.0d+00)

        do ig = 1,ng
          do ir = 1,nr
            do it = 1,nt
              do il = 1,nul
                ff(il,it,ir,ig) = zero
              enddo
              fn(it,ir,ig) = zero
            enddo
          enddo
        enddo
        ig = ng9
          do ir = 1,nr9
            do it = 1,nt
              do il = 1,nul
                ff(il,it,ir,ig) = zero
              enddo
              fn(it,ir,ig) = zero
            enddo
          enddo

        return
      end


      subroutine zerog(nul)
!	zero gg and gn

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mu = 64)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg), &
                ff(mu,mt,mr,mg),fn(mt,mr,mg)
        parameter (zero = 0.0d+00)

        do ig = 1,ng
          do ir = 1,nr
            do it = 1,nt
              do il = 1,nul
                gg(il,it,ir,ig) = zero
              enddo
              gn(it,ir,ig) = zero
            enddo
          enddo
        enddo
        ig = ng9
          do ir = 1,nr9
            do it = 1,nt
              do il = 1,nul
                gg(il,it,ir,ig) = zero
              enddo
              gn(it,ir,ig) = zero
            enddo
          enddo

        return
      end



      subroutine copyarr(arr)
!	copy array between grids

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        common /mpis/mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        dimension arr(mt,mr,mg)
        logical badnum,doprint
        parameter (zero = 0.0d+00)
        parameter (two = 2.0d+00)

        doprint = (myrank.eq.0)
        nr2 = nr/2
!	copy from expanded grid onto grid 1
        ig = 1
        jg = ng9
        xn9 = nn9**2
        nn8 = nn9-1
        nt8 = nt-nt/nn9+1
        do ir = 1,nr
          jr2 = nn9*ir
          jr1 = jr2-nn8
          do it = 1,nt,nn9
            suma = zero
            do jr = jr1,jr2
              do jt = it,it+nn8
                suma = suma+arr(jt,jr,jg)
              enddo
            enddo
            if(badnum(suma)) then
              if(doprint) print '(e10.2," in copyarr at",3i6)', &
                suma,it,jr1,jg
              suma = zero
            endif
            arr(it/nn9+nt8,ir,ig) = suma/xn9
          enddo
        enddo
!	copy to inner half of next grid
        do ig = 2,ng
          jg = ig-1
          do ir = 1,nr2
            jr2 = 2*ir
            jr1 = jr2-1
            do it = 1,nt
              arr(it,ir,ig) = (arr(it,jr1,jg)+arr(it,jr2,jg))/two
              if(badnum(arr(it,ir,ig))) then
                if(doprint) print '(e10.2," in copyarr at",3i6)', &
                  arr(it,ir,ig),it,jr1,jg
                arr(it,ir,ig) = zero
              endif
            enddo
          enddo
        enddo

        return
      end



      subroutine copytau(arr)
!	copy array between grids without averaging

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        common /mpis/mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        dimension arr(mt,mr,mg)
        logical badnum,doprint
        parameter (zero = 0.0d+00)

        doprint = (myrank.eq.0)
        nr2 = nr/2
!	copy from expanded grid onto grid 1
        ig = 1
        jg = ng9
        xn9 = nn9
        nn8 = nn9-1
        nt8 = nt-nt/nn9+1
        do ir = 1,nr
          jr = nn9*ir
          do it = 1,nt,nn9
            suma = zero
            do jt = it,it+nn8
              suma = suma+arr(jt,jr,jg)
            enddo
            if(badnum(suma)) then
              if(doprint) print '(e10.2," in copytau at",3i6)', &
                suma,it,jr,jg
              suma = zero
            endif
            arr(it/nn9+nt8,ir,ig) = suma/xn9
          enddo
        enddo
!	copy to inner half of next grid
        do ig = 2,ng
          jg = ig-1
          do ir = 1,nr2
            jr = 2*ir
            do it = 1,nt
              arr(it,ir,ig) = arr(it,jr,jg)
            enddo
          enddo
        enddo

        return
      end



      subroutine copywta(arr,wtarr)
!	copy between grids with weighting

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (zero = 0.0d+00)
        parameter (two = 2.0d+00)
        common /mpis/mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        dimension arr(mt,mr,mg),wtarr(mt,mr,mg)
        logical badnum,doprint

        doprint = (myrank.eq.0)
        nr2 = nr/2
        nn8 = nn9-1
        nt8 = nt-nt/nn9+1
        xn9 = nn9**2
!	copy from expanded grid onto grid 1
        ig = 1
        jg = ng9
        do ir = 1,nr
          jr2 = nn9*ir
          jr1 = jr2-nn8
          do it = 1,nt,nn9
            suma = zero
            sumwt = zero
            do jr = jr1,jr2
              do jt = it,it+nn8
                suma = suma+arr(jt,jr,jg)*wtarr(jt,jr,jg)
                sumwt = sumwt+wtarr(jt,jr,jg)
              enddo
              if(badnum(suma)) then
                if(doprint) print '(e10.2," in copywta at",3i6)', &
                  suma,it,jr1,jg
                suma = zero
              endif
            enddo
            if(sumwt.gt.zero) then
              arr(it/nn9+nt8,ir,ig) = suma/sumwt
            else
              suma = zero
              do jr = jr1,jr2
                do jt = it,it+nn8
                  suma = suma+arr(jt,jr,jg)
                enddo
                if(badnum(suma)) then
                  if(doprint) print '(e10.2," in copywta at",3i6)', &
                    suma,it,jr1,jg
                  suma = zero
                endif
              enddo
              arr(it/nn9+nt8,ir,ig) = suma/xn9
            endif
          enddo
        enddo
!	copy to inner half of next grid
        do ig = 2,ng
          jg = ig-1
          do ir = 1,nr2
            jr2 = 2*ir
            jr1 = jr2-1
            do it = 1,nt
              sumwt = wtarr(it,jr1,jg)+wtarr(it,jr2,jg)
              if(sumwt.gt.zero) then
                arr(it,ir,ig) = (arr(it,jr1,jg)*wtarr(it,jr1,jg) &
                  +arr(it,jr2,jg)*wtarr(it,jr2,jg))/sumwt
              else
                arr(it,ir,ig) = (arr(it,jr1,jg)+arr(it,jr2,jg))/two
              endif
              if(badnum(arr(it,ir,ig))) then
                if(doprint) print '(e10.2," in copywta at",3i6)', &
                  arr(it,ir,ig),it,jr1,jg
                arr(it,ir,ig) = zero
              endif
            enddo
          enddo
        enddo

        return
      end



      subroutine copyarrr(arrr,mul,nul)
!	copy 4-d array between grids

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        common /mpis/mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        logical badnum,doprint
        dimension arrr(mul,mt,mr,mg)
        parameter (zero = 0.0d+00)
        parameter (two = 2.0d+00)

        doprint = (myrank.eq.0)
        nr2 = nr/2
        nn8 = nn9-1
        xn9 = nn9**2
        nt8 = nt-nt/nn9+1
!	copy from expanded grid onto grid 1
        ig = 1
        jg = ng9
        do ir = 1,nr
          jr2 = nn9*ir
          jr1 = jr2-nn8
          do it = 1,nt,nn9
            do il = 1,nul
              suma = zero
              do jr = jr1,jr2
                do jt = it,it+nn8
                  suma = suma+arrr(il,jt,jr,jg)
                enddo
              enddo
              if(badnum(suma)) then
                if(doprint) print '(e10.2," in copyarrr at",4i6)', &
                  suma,il,it,jr1,jg
                suma = zero
              endif
              arrr(il,it/nn9+nt8,ir,ig) = suma/xn9
            enddo
          enddo
        enddo
!	copy to inner half of next grid
        do ig = 2,ng
          jg = ig-1
          do ir = 1,nr2
            jr2 = 2*ir
            jr1 = jr2-1
            do it = 1,nt
              do il = 1,nul
                arrr(il,it,ir,ig) &
                  = (arrr(il,it,jr1,jg)+arrr(il,it,jr2,jg))/two
                if(badnum(arrr(il,it,ir,ig))) then
                  if(doprint) print '(e10.2," in copyarrr at",4i6)', &
                    arrr(il,it,ir,ig),il,it,jr1,jg
                  arrr(il,it,ir,ig) = zero
                endif
              enddo
            enddo
          enddo
        enddo

        return
      end



      subroutine plotcon(con,logplot,monfile)
!	make r-th contour plots

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (zero = 0.0d+00)
        parameter (two = 2.0d+00)
        parameter (four = 4.0d+00)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        dimension con(mt,mr,mg)
        logical logplot,monplot
        character*16 lockline,monfile
        dimension conmap(2*mt)

        monplot = (index(monfile,' ').gt.1)
        nr2 = nr/2

   86   open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=90)
        read(8,'(a16)') lockline
        close(unit=8)
        if(index(lockline,'reading').gt.0) then
          call sleep(1)
          goto 86
        endif
        open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=90)
        write(8,'("writing")')
        close(unit=8)
        open(unit=9,file='/home/r900-1/lacy/plots/contours', &
          access='direct',err=90,recl=2*nt)
        if(monplot) open(unit=10,file=monfile,access='direct',recl=2*nt)
 
!	plot map of grid ng9
        irec = 0
        ig = ng9
        if(nr9.le.nt) then
          nxm = nt
          nym = nr9
          do ir = 1,nr9
            imap = 0
            do it = 1,nt
              imap = imap+1
              if(logplot) then
                if(con(it,ir,ig).gt.zero) then
                  conmap(imap) = slog10(con(it,ir,ig))
                else
                  conmap(imap) = zero
                endif
              else
                conmap(imap) = con(it,ir,ig)
              endif
            enddo  !  it
            irec = irec+1
            write(unit=9,rec=irec) (conmap(it),it=1,nt)
          enddo  !  ir
        elseif(nr9.le.(2*nt)) then
          nxm = 2*nt
          nym = nr9
          do ir = 1,nr9
            imap = 0
            do it = 1,nt
              imap = imap+2
              if(logplot) then
                if(con(it,ir,ig).gt.zero) then
                  conmap(imap) = slog10(con(it,ir,ig))
                else
                  conmap(imap) = zero
                endif
              else
                conmap(imap) = con(it,ir,ig)
              endif
              if(imap.gt.2) then
                conmap(imap-1) = (conmap(imap-2)+conmap(imap))/two
              else
                conmap(1) = conmap(2)
              endif
            enddo  !  it
            irec = irec+1
            write(unit=9,rec=irec) (conmap(it),it=1,nt)
            irec = irec+1
            write(unit=9,rec=irec) (conmap(it),it=nt+1,2*nt)
            if(monplot) then
              write(unit=10,rec=irec-1) (conmap(it),it=1,nt)
              write(unit=10,rec=irec) (conmap(it),it=nt+1,2*nt)
            endif
          enddo  !  ir
        else
          nxm = 2*nt
          nym = nr9/2
          do ir = 1,nr9,2
            imap = 0
            do it = 1,nt
              imap = imap+2
              if(logplot) then
                prod = con(it,ir,ig)*con(it,ir+1,ig)
                if(prod.gt.zero) then
                  conmap(imap) = slog10(prod)/two
                else
                  conmap(imap) = zero
                endif
              else
                conmap(imap) = (con(it,ir,ig)+con(it,ir+1,ig))/two
              endif
              if(imap.gt.2) then
                conmap(imap-1) = (conmap(imap-2)+conmap(imap))/two
              else
                conmap(1) = conmap(2)
              endif
            enddo  !  it
            irec = irec+1
            write(unit=9,rec=irec) (conmap(it),it=1,nt)
            irec = irec+1
            write(unit=9,rec=irec) (conmap(it),it=nt+1,2*nt)
            if(monplot) then
              write(unit=10,rec=irec-1) (conmap(it),it=1,nt)
              write(unit=9,rec=irec) (conmap(it),it=nt+1,2*nt)
            endif
          enddo  !  ir
        endif

        close(unit=9)
        open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=90)
        write(8,'("contours")')
        write(8,'(3i8)') nxm,nym,1
        close(unit=8)
        call waitasec(.false.,4)

   88   open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=90)
        read(8,'(a16)') lockline
        close(unit=8)
        if(index(lockline,'reading').gt.0) then
          call sleep(1)
          goto 88
        endif
        open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=90)
        write(8,'("writing")')
        close(unit=8)
        open(unit=9,file='/home/r900-1/lacy/plots/contours', &
          access='direct',err=90,recl=2*nt)
 
!	plot map of grids 1-ng
        jrec = irec
        irec = 0
        if((ng*nr/2).le.(2*nt)) then
          nxm = 2*nt
          nym = ng*nr2
          do ig = 1,ng
            do ir = nr2+1,nr
              imap = 0
              do it = 1,nt
                imap = imap+2
                if(logplot) then
                  if(con(it,ir,ig).gt.zero) then
                    conmap(imap) = slog10(con(it,ir,ig))
                  else
                    conmap(imap) = zero
                  endif
                else
                  conmap(imap) = con(it,ir,ig)
                endif
                if(imap.gt.2) then
                  conmap(imap-1) = (conmap(imap-2)+conmap(imap))/two
                else
                  conmap(1) = conmap(2)
                endif
              enddo  !  it
              irec = irec+1
              write(unit=9,rec=irec) (conmap(it),it=1,nt)
              irec = irec+1
              write(unit=9,rec=irec) (conmap(it),it=nt+1,2*nt)
              if(monplot) then
                jrec = jrec+1
                write(unit=10,rec=jrec) (conmap(it),it=1,nt)
                jrec = jrec+1
                write(unit=10,rec=jrec) (conmap(it),it=nt+1,2*nt)
              endif
            enddo  !  ir
          enddo  !  ig
        elseif((ng*nr/2).le.(4*nt)) then
          nxm = 2*nt
          nym = ng*nr2/2
          do ig = 1,ng
            do ir = nr2+1,nr,2
              imap = 0
              do it = 1,nt
                imap = imap+2
                if(logplot) then
                  prod = con(it,ir,ig)*con(it,ir+1,ig)
                  if(prod.gt.zero) then
                    conmap(imap) = slog10(prod)/two
                  else
                    conmap(imap) = zero
                  endif
                else
                  conmap(imap) = (con(it,ir,ig)+con(it,ir+1,ig))/two
                endif
                if(imap.gt.2) then
                  conmap(imap-1) = (conmap(imap-2)+conmap(imap))/two
                else
                  conmap(1) = conmap(2)
                endif
              enddo  !  it
              irec = irec+1
              write(unit=9,rec=irec) (conmap(it),it=1,nt)
              irec = irec+1
              write(unit=9,rec=irec) (conmap(it),it=nt+1,2*nt)
              if(monplot) then
                jrec = jrec+1
                write(unit=10,rec=jrec) (conmap(it),it=1,nt)
                jrec = jrec+1
                write(unit=10,rec=jrec) (conmap(it),it=nt+1,2*nt)
              endif
            enddo  !  ir
          enddo  !  ig
        else
          nxm = 2*nt
          nym = ng*nr2/4
          do ig = 1,ng
            do ir = nr2+1,nr,4
              imap = 0
              do it = 1,nt
                imap = imap+2
                if(logplot) then
                  prod = con(it,ir,ig)*con(it,ir+1,ig) &
                      *con(it,ir+2,ig)*con(it,ir+3,ig)
                  if(prod.gt.zero) then
                    conmap(imap) = slog10(prod)/four
                  else
                    conmap(imap) = zero
                  endif
                else
                  conmap(imap) = (con(it,ir,ig)+con(it,ir+1,ig) &
                    +con(it,ir+2,ig)+con(it,ir+3,ig))/four
                endif
                if(imap.gt.2) then
                  conmap(imap-1) = (conmap(imap-2)+conmap(imap))/two
                else
                  conmap(1) = conmap(2)
                endif
              enddo  !  it
              irec = irec+1
              write(unit=9,rec=irec) (conmap(it),it=1,nt)
              irec = irec+1
              write(unit=9,rec=irec) (conmap(it),it=nt+1,2*nt)
              if(monplot) then
                jrec = jrec+1
                write(unit=10,rec=jrec) (conmap(it),it=1,nt)
                jrec = jrec+1
                write(unit=10,rec=jrec) (conmap(it),it=nt+1,2*nt)
              endif
            enddo  !  ir
          enddo  !  ig
        endif

        close(unit=9)
        if(monplot) close(unit=10)
        open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=90)
        write(8,'("contours")')
        write(8,'(3i8)') nxm,nym,1
        close(unit=8)

 90     return
      end



      subroutine plotim(frame,mx,nx,monfile)
!	make greyscale plot of an image frame

        implicit real*8 (a-h,o-z)
        dimension frame(mx,mx)
        character*16 lockline,monfile
        logical monplot

        monplot = (index(monfile,' ').gt.1)
   88   open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=90)
        read(8,'(a16)') lockline
        close(unit=8)
        if(index(lockline,'reading').gt.0) then
          call sleep(1)
          goto 88
        endif
        open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=90)
        write(8,'("writing")')
        close(unit=8)
        open(unit=9,file='/home/r900-1/lacy/plots/images', &
          access='direct',recl=2*nx)
        if(monplot) open(unit=10,file=monfile,access='direct',recl=2*nx)

        do iy = 1,nx
          write(unit=9,rec=iy) (frame(ix,iy),ix=1,nx)
          if(monplot) write(unit=10,rec=iy) (frame(ix,iy),ix=1,nx)
        enddo

        close(unit=9)
        if(monplot) close(unit=10)
        open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=90)
        write(8,'("images")')
        write(8,'(3i8)') nx,nx,1
        close(unit=8)

   90   return
      end


      subroutine plotggfs(gg,jg,ju)

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mu = 64)
        parameter (ma = 32)
        parameter (mc = 1024)
        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (hun = 1.0d+02)

        common /mpis/mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu), &
                phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg), &
                opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg), &
                rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        dimension gg(mu,mt,mr,mg)
        dimension frame(mc,mc)
        logical badnum

        nc = mc
        alog2 = log(two)
        call startgr(jg,ig,kr,lr,lt,dr,dcost)
        if(ig.gt.ng) then
!          dxm = (nr9*dr9)/(two*nc)
!          x0 = (nc-half)*dxm
          dxm = (nr9*dr9)/nc
          x0 = -half*dxm
          zstr = nn9
        else
          dxm = (nr*dr)/nc
          x0 = -half*dxm
          zstr = one
        endif
        rnr = nr*dr0
        summap = zero
        summi = zero
        do ix = 1,nc
          do iz = 1,nc
            x = x0+ix*dxm
            y = zero
            z = (iz-nc/2-half)*dxm
!	stretch z scale
            z = z/zstr
            r = sqrt(x**2+y**2+z**2)
            costh = abs(z/r)
            if(ig.gt.ng) then
              it = nt-nint(costh/dcost-half)
              if((it.ge.1).and.(it.le.nt)) then
                call startgr(jg,kg,kr,lr,lt,drgr,dcosgr)
                ir = r/dr9+one
                ir = min(nr9,max(1,ir))
              else
                igr = slog(r/rnr)/alog2+one
                igr = min(ng,max(1,igr))
                call startgr(igr,kg,kr,lr,lt,drgr,dcosgr)
                ir = r/drgr+one
                it = nt-nint(costh/dcosgr-half)
                if((it.lt.1).or.(it.gt.nt).or.(ir.lt.1).or.(ir.gt.nr)) &
                  then
                  frame(ix,iz) = one
                  goto 240
                endif
              endif
            else
              igr = slog(r/rnr)/alog2+one
              igr = min(ng,max(1,igr))
              call startgr(igr,kg,kr,lr,lt,drgr,dcosgr)
              ir = r/drgr+one
              it = nt-nint(costh/dcosgr-half)
              if((it.lt.1).or.(it.gt.nt).or.(ir.lt.1).or.(ir.gt.nr)) &
                then
                frame(ix,iz) = one
                goto 240
              endif
            endif
            if(r.le.rstar) then
              frame(ix,iz) = zero
              goto 240
            elseif(r.lt.(hun*rstar)) then
              ostar = (one-sqrt(one-(rstar/r)**2))/two
            else
              ostar = (rstar/(two*r))**2
            endif
            if(ju.eq.0) then
              ggsum = zero
              fssum = zero
              do iu = 1,nu
                ggsum = ggsum+opa(iu,nt,nr,1)*gg(iu,it,ir,kg)*dud(iu)
                fssum = fssum+opa(iu,nt,nr,1)*bus(iu)*dud(iu)
              enddo
              if(fssum.gt.zero) then
                ggfs = ggsum/(fssum*ostar)
              else
                ggfs = one
              endif
            elseif(ju.lt.0) then
              ia = abs(ju)
              ggsum = zero
              fssum = zero
              do iu = 1,nu
                ggsum = ggsum+phab(iu,ia)*gg(iu,it,ir,kg)*dud(iu)
                fssum = fssum+phab(iu,ia)*bus(iu)*dud(iu)
              enddo
              if(fssum.gt.zero) then
                ggfs = ggsum/(fssum*ostar)
              else
                ggfs = one
              endif
            else
              ggfs = gg(ju,it,ir,kg)/(bus(ju)*ostar)
            endif
            if(badnum(ggfs)) then
              print '(" in ggfs:",5i6,1p,2e10.2)', &
                ix,iz,it,ir,kg,gg(ju,it,ir,kg),bus(iu)
              ggfs = zero
            endif
            summap = summap+ggfs
            summi = summi+one
            frame(ix,iz) = sqrt(ggfs)
  240     enddo
        enddo

        if(ju.lt.0) then
          print'(" sqrt(gg(ia=",i2,")/stellar) with z stretch =",f4.1)', &
          ia,zstr
        elseif(ju.eq.0) then
          print'(" sqrt(gg(condensed)/stellar) with z stretch =",f4.1)', &
          zstr
        else
          print'(" sqrt(gg(iu=",i2,")/stellar) with z stretch =",f4.1)', &
          ju,zstr
        endif
!        print '(" mean gg/fs =",1p,e10.3)', summap/summi
        call plotim(frame,mc,nc,'')

        return
      end


      subroutine plotxyim(pl,jg)

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mc = 1024)
        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)

        common /mpis/mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        dimension pl(mt,mr,mg)
        dimension frame(mc,mc)

        nc = mc
        alog2 = log(two)
        call startgr(jg,ig,kr,lr,lt,dr,dcost)
        if(ig.gt.ng) then
          rnr = kr*dr
          dxm = rnr/(2*nc)
          x0 = (nc-half)*dxm
        else
          dxm = dr
          rnr = nr*dr0
          x0 = -half*dxm
        endif
        summap = zero
        summi = zero
        do ix = 1,nc
          do iz = 1,nc
            x = x0+ix*dxm
            y = zero
            z = (iz-nc/2-half)*dxm
            r = sqrt(x**2+y**2+z**2)
            costh = abs(z/r)
            if(ig.gt.ng) then
              it = nt-nint(costh/dcost-half)
              if((it.ge.1).and.(it.le.nt)) then
                igr = ig
                drgr = dr9
                ir = r/dr9+one
                ir = min(nr9,max(1,ir))
              else
                igr = slog(r/rnr)/alog2+one
                igr = min(ng,max(1,igr))
                call startgr(igr,kg,kr,lr,lt,dr,dcosgr)
                drgr = dr
                ir = r/drgr+one
                ir = min(nr,max(1,ir))
                it = nt-nint(costh/dcosgr-half)
                if((it.lt.1).or.(it.gt.nt)) then
                  frame(ix,iy) = zero
                  goto 240
                endif
              endif
            else
              igr = slog(r/rnr)/alog2+one
              igr = min(ng,max(1,igr))
              call startgr(igr,kg,kr,lr,lt,dr,dcosgr)
              drgr = dr
              ir = r/drgr+one
              ir = min(nr,max(1,ir))
              it = nt-nint(costh/dcosgr-half)
              if((it.lt.1).or.(it.gt.nt)) then
                frame(ix,iy) = zero
                goto 240
              endif
            endif
            if(r.le.rstar) then
              frame(ix,iz) = zero
              goto 240
            elseif(r.lt.(100.*rstar)) then
              ostar = (one-sqrt(one-(rstar/r)**2))/2.
            else
              ostar = 0.25*(rstar/r)**2
            endif
!            frame(ix,iz) = pl(it,ir,kg)/ostar
            frame(ix,iz) = pl(it,ir,kg)*(r/rstar)
            summap = summap+frame(ix,iz)
            summi = summi+one
            frame(ix,iz) = max(frame(ix,iz),zero)
  240     enddo
        enddo

        call plotim(frame,mc,nc,'')

        return
      end


      subroutine plotspec(spec,mw,nw,nspec,irot,dw,monfile)
!	plot spectrum

        implicit real*8 (a-h,o-z)
        dimension spec(mw,nspec)
        character*16 lockline,monfile
        logical monplot

        monplot = (index(monfile,' ').gt.1)
   88   open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=90)
        read(8,'(a16)') lockline
        close(unit=8)
        if(index(lockline,'reading').gt.0) then
          call sleep(1)
          goto 88
        endif
        open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=90)
        write(8,'("writing")')
        close(unit=8)
        open(unit=9,file='/home/r900-1/lacy/plots/specs', &
          access='direct',recl=2*nw)
        if(monplot) open(unit=10,file=monfile,access='direct',recl=2*nw)

        do ispec = 1,nspec
          write(unit=9,rec=ispec) (spec(iw,ispec),iw=1,nw)
          jspec = (irot-1)*nspec+ispec
          if(monplot) write(unit=10,rec=jspec) (spec(iw,ispec),iw=1,nw)
        enddo

        close(unit=9)
        if(monplot) close(unit=10)
        open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
          access='sequential',err=90)
        write(8,'("specs")')
        write(8,'(2i8,1p,e12.3)') nw,nspec,dw
        close(unit=8)

   90   return
      end



      function bw(w,t)
!	black body function

        implicit real*8 (a-h,o-z)
        parameter (hc2 = 5.955d-06)
        parameter (hck = 1.4388d+00)
        parameter (big = 5.0d+02)
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)

        arg = hck*w/t
        if(arg.le.zero) then
          bw = zero
        elseif(arg.lt.big) then
          bw = two*hc2*w**3/(sexp(arg)-one)
        else
          bw = zero
        endif

        return
      end


      function texp(arg)
!	exp(arg) using look-up table

        implicit real*8 (a-h,o-z)
        common /exptab/ expi(10239),too
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (hun = 1.0d+02)
        parameter (big = 1.0d+04)

        harg = abs(hun*arg)
        if(harg.lt.big) then
          iarg = int(harg)
          darg = (harg-iarg)/hun
          texp = (one-darg+darg**2/too)*expi(iarg+1)
          if(arg.gt.zero) texp = one/texp
        elseif(arg.lt.zero) then
          texp = zero
        else
          texp = one/expi(10000)
        endif

        return
 
      entry expinit

        do iarg = 1,10000
          expi(iarg) = sexp(-(iarg-1)/hun)
        enddo
        darg = one/hun
        too = darg**2/(exp(-darg)-one+darg)
        expinit = one

        return
      end


      function sexp(arg)
!	safe and efficient exp function

        implicit real*8 (a-h,o-z)
        parameter (big = 4.6d+02)
        parameter (small = 1.0d-03)
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        save nan
        if((nan.lt.0).or.(nan.gt.8)) nan = 0

        if(arg.lt.(-big)) then
          sexp = zero
        elseif(arg.gt.big) then
          sexp = 1.0d+200
        elseif(abs(arg).lt.small) then
          sexp = 1.0+arg+arg**2/two
        elseif((arg.lt.zero).or.(arg.gt.zero)) then
          sexp = exp(arg)
        else  !  NaN
          print '(e10.2," passed to sexp",i6)', arg,nan
          nan = nan+1
          if(nan.gt.8) stop
          sexp = one
        endif

        return
      end


      function slog(arg)
!	safe log function

        implicit real*8 (a-h,o-z)
        parameter (tiny = 1.0d-200)
        parameter (big = 1.0d+200)
        parameter (small = 1.0d-03)
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        save nan
        if((nan.lt.0).or.(nan.gt.8)) nan = 0

        if(arg.le.zero) then
          slog = -500.
        elseif(arg.lt.tiny) then
          slog = -460.
        elseif(arg.gt.big) then
          slog = 460.
        elseif(abs(arg-one).lt.small) then
          slog = arg-one
          slog = slog-slog**2/two
        elseif(arg.gt.tiny) then
          slog = log(arg)
        else  !  NaN
          if(nan.le.8) print '(e10.2," passed to slog",i6)', arg
          nan = nan+1
          slog = zero
        endif

        return
      end


      function slog10(arg)
!	safe log_10 function

        implicit real*8 (a-h,o-z)
        parameter (tiny = 1.0d-200)
        parameter (big = 1.0d+200)
        parameter (small = 1.0d-03)
        parameter (aloge = 0.4342945d+00)
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        save nan
        if((nan.lt.0).or.(nan.gt.8)) nan = 0

        if(arg.le.zero) then
          slog10 = -250.
        elseif(arg.lt.tiny) then
          slog10 = -200.
        elseif(arg.gt.big) then
          slog10 = 200.
        elseif(abs(arg-one).lt.small) then
          slog = arg-one
          slog10 = aloge*(slog-slog**2/two)
        elseif(arg.gt.tiny) then
          slog10 = log10(arg)
        else  !  NaN
          if(nan.le.8) print '(e10.2," passed to slog10",i6)', arg
          nan = nan+1
          slog10 = zero
        endif

        return
      end


      subroutine waitasec(wait,iwait)

        implicit real*8 (a-h,o-z)
        common /mpis/mysize,myrank
        logical wait
        character*16 lockline

        if(wait) then
          call sleep(iwait-1)
  100     call sleep(1)
          open(unit=8,file='/home/r900-1/lacy/plots/plotlock', &
            access='sequential',err=90)
          read(8,'(a16)') lockline
          close(unit=8)
          if(index(lockline,'empty').eq.0) goto 100
        else
          call sleep(iwait)
        endif

   90   return
      end


      subroutine savemodel

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        parameter (mw = 1024)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu), &
                phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg), &
                opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg), &
                rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml), &
                wl(ml),dl(ml,mt,mr,mg),sl(ml,mt,mr,mg),opl(ml,mt,mr,mg)
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg), &
                ff(mu,mt,mr,mg),fn(mt,mr,mg)
        parameter (half = 0.5d+00)
        parameter (pi = 3.14159265d+00)

        open(unit=9,file='diskmodel',access='direct',recl=2*nt)
        irec = 0
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            irec = irec+1
            write(unit=9,rec=irec) (dg(it,ir,ig),it=1,nt)
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            irec = irec+1
            write(unit=9,rec=irec) (tg(it,ir,ig),it=1,nt)
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            irec = irec+1
            write(unit=9,rec=irec) (td(it,ir,ig),it=1,nt)
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            irec = irec+1
            write(unit=9,rec=irec) (opl(ml,it,ir,ig),it=1,nt)
          enddo
        enddo
        close(unit=9)
        print '(" wrote to diskmodel",2i8," byte records")', irec,8*nt
        irec = 0
        open(unit=9,file='dustmodel',access='direct',recl=2*nu)
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            do it = 1,nt
              irec = irec+1
              write(unit=9,rec=irec) (opa(iu,it,ir,ig),iu=1,nu)
            enddo
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            do it = 1,nt
              irec = irec+1
              write(unit=9,rec=irec) (swd(iu,it,ir,ig),iu=1,nu)
            enddo
          enddo
        enddo
        close(unit=9)
        print '(" wrote to dustmodel",2i8," byte records")', irec,8*nu

        open(unit=10,file='disk.out')
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            r = dr*(ir-half)
            do it = 1,nt
              z = r*dcosth*(nt-it+half)
              sumd = zero
              ddi = dd(it,ir,ig)
              dgi = dg(it,ir,ig)
              do ia = 1,na
                sumd = sumd+condens(ia,it,ir,ig)*sexp(-ad(ia)*ddi) &
                       *dgi*popd(ia)*pi*ad(ia)**2
              enddo
              write(10,'(3i6,1p,7e10.3)') igr,ir,it,r,z, &
                dg(it,ir,ig),tg(it,ir,ig),td(it,ir,ig),gguv(it,ir,ig), &
                sumd
            enddo
          enddo
        enddo

        return
      end


      subroutine readmodel

        implicit real*8 (a-h,o-z)
#include "mpif.h"
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        parameter (mw = 1024)
        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu), &
                phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg), &
                opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg), &
                rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml), &
                wl(ml),dl(ml,mt,mr,mg),sl(ml,mt,mr,mg),opl(ml,mt,mr,mg)
        logical badnum,doprint

      if(myrank.eq.0) then
        irec = 0
        open(unit=9,file='diskmodel',access='direct',recl=2*nt,err=98)
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            irec = irec+1
            read(unit=9,rec=irec,err=98) (dg(it,ir,ig),it=1,nt)
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            irec = irec+1
            read(unit=9,rec=irec,err=98) (tg(it,ir,ig),it=1,nt)
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            irec = irec+1
            read(unit=9,rec=irec,err=98) (td(it,ir,ig),it=1,nt)
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            irec = irec+1
            read(unit=9,rec=irec,err=98) (opl(ml,it,ir,ig),it=1,nt)
          enddo
        enddo
        close(unit=9)

        irec = 0
        open(unit=9,file='dustmodel',access='direct',recl=2*nu,err=99)
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            do it = 1,nt
              irec = irec+1
              read(unit=9,rec=irec,err=99) (opa(iu,it,ir,ig),iu=1,nu)
            enddo
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            do it = 1,nt
              irec = irec+1
              read(unit=9,rec=irec,err=99) (swd(iu,it,ir,ig),iu=1,nu)
            enddo
          enddo
        enddo
        close(unit=9)

        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            do it = 1,nt
              if(doprint.and.badnum(dg(it,ir,ig))) &
                print '("dg(",3i4,") =",e10.2)', it,ir,ig,dg(it,ir,ig)
              if(doprint.and.badnum(tg(it,ir,ig))) &
                print '("tg(",3i4,") =",e10.2)', it,ir,ig,tg(it,ir,ig)
              if(doprint.and.badnum(td(it,ir,ig))) &
                print '("td(",3i4,") =",e10.2)', it,ir,ig,td(it,ir,ig)
              if(doprint.and.badnum(opl(ml,it,ir,ig))) print &
                '("opl(",3i4,") =",e10.2)', it,ir,ig,opl(ml,it,ir,ig)
              do iu = 1,nu
                if(doprint.and.badnum(opa(iu,it,ir,ig))) &
                  print '("opa(",4i4,") =",e10.2)', &
                  iu,it,ir,ig,opa(iu,it,ir,ig)
                if(doprint.and.badnum(swd(iu,it,ir,ig))) &
                  print '("swd(",4i4,") =",e10.2)', &
                  iu,it,ir,ig,swd(iu,it,ir,ig)
              enddo
            enddo
          enddo
        enddo
      endif

      if(mysize.gt.0) then
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            call mpi_bcast(dg(1,ir,ig),nt, &
              MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
            call mpi_bcast(tg(1,ir,ig),nt, &
              MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
            call mpi_bcast(td(1,ir,ig),nt, &
              MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
            do it = 1,nt
              call mpi_bcast(opl(ml,it,ir,ig),1, &
                MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
              call mpi_bcast(opa(1,it,ir,ig),nu, &
                MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
              call mpi_bcast(swd(1,it,ir,ig),nu, &
                MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
            enddo
          enddo
        enddo
      endif

        call copyarr(dg)
        call copyarr(tg)
        call copyarr(td)
        call copyarrr(opl,ml,ml)
        call copyarrr(opa,mu,nu)
        call copyarrr(swd,mu,nu)

        return

   98   print '(" Error reading diskmodel",5i8)', ir,ig,irec,8*nt
        stop

   99   print '(" Error reading dustmodel",5i8)', it,ir,ig,irec,8*nu
        stop

      end



      subroutine readparms

        implicit real*8 (a-h,o-z)
        parameter (mc = 256)

        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /heat/ gmmdot,acom,pxr,axr,bxr,epspe,agas,sigma, &
                colin,colup,xco,xwat
        common /angls/ dcost,dcos9,rot,irot
        common /parmz/ rgrid1,r0,r1,sigma0,sigma1,ocol,pcol,qcol, &
                avisc,xmstar,xlstar,tuv,xuv,tsubl,gsubl,tburn,xburn, &
                coswind,xc,xsi,ac,asi,enhanc,xmdot,wmdot,vwind, &
                xco0,xwat0,tfrco,tfrwat,xmol,tfrmol,tdiss,coldiss, &
                kit,lit,mit,nnt,lx,lrot,reuse,lte
        common /cube/ cube(mc,mc,mc),ngc
        logical doprint,reuse,lte
        character*80 line,parn,parv

        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (half = 0.5d+00)
        parameter (quart = 0.25d+00)
        parameter (pi = 3.14159265d+00)
        parameter (au = 1.496d+13)
        parameter (rsun = 6.96d+10)
        parameter (xlsun = 3.90d+33)
        parameter (xmsun = 1.99d+33)
!	mass / gas particle
        parameter (xmg = 4.0d-24)
        parameter (vkms = 1.0d+05)

        doprint = (myrank.eq.0)
        open(unit=12,file='frisco.par')

  100   read(unit=12,fmt='(a72)',err=90,end=90) line
        if(index(line,'c').eq.1) then
          continue
        elseif(index(line,'=').gt.0) then
          if(doprint) print '(a72)', line(1:72)
          npar = index(line,'=')-1
          parn = line(1:npar)
          npar = npar+2
          parv = line(npar:72)
          if(index(parn,'reuse').gt.1) then
            reuse = (index(parv,'t').gt.1)
          elseif(index(parn,'lte').gt.1) then
            lte = (index(parv,'t').gt.1)
          elseif(index(parn,'rgrid1').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) rgrid1
            rgrid1 = rgrid1*au
          elseif(index(parn,'r0 ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) r0
            r0 = r0*au
          elseif(index(parn,'r1 ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) r1
            r1 = r1*au
          elseif(index(parn,'sigma0').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) sigma0
          elseif(index(parn,'sigma1').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) sigma1
          elseif(index(parn,'pcol').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) pcol
          elseif(index(parn,'qcol').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) qcol
          elseif(index(parn,'ocol').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) ocol
          elseif(index(parn,'avisc').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) avisc
          elseif(index(parn,'rstar').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) rstar
            rstar = rstar*rsun
          elseif(index(parn,'mstar').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xmstar
          elseif(index(parn,'lstar').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xlstar
          elseif(index(parn,'tstar').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tstar
          elseif(index(parn,'tuv ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tuv
          elseif(index(parn,'xuv ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xuv
          elseif(index(parn,'tsubl').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tsubl
          elseif(index(parn,'gsubl').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) gsubl
          elseif(index(parn,'tburn').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tburn
          elseif(index(parn,'xburn').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xburn
          elseif(index(parn,'enhanc').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) enhanc
          elseif(index(parn,'wmdot').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) wmdot
          elseif(index(parn,'mdot').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xmdot
          elseif(index(parn,'acom').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) acom
          elseif(index(parn,'pxr ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) pxr
          elseif(index(parn,'axr ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) axr
          elseif(index(parn,'bxr ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) bxr
          elseif(index(parn,'epspe').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) epspe
          elseif(index(parn,'xco').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xco0
          elseif(index(parn,'xwat').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xwat0
          elseif(index(parn,'tfrco').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tfrco
          elseif(index(parn,'tfrwat').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tfrwat
          elseif(index(parn,'coswind').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) coswind
          elseif(index(parn,'vwind').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) vwind
          elseif(index(parn,'dvturb0').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) dvturb0
            dvturb0 = dvturb0*vkms
          elseif(index(parn,'dvturb1').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) dvturb1
            dvturb1 = dvturb1*vkms
          elseif(index(parn,'dv ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) dv
            dv = dv*vkms
          elseif(index(parn,'xmol').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xmol
          elseif(index(parn,'coldiss').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) coldiss
          elseif(index(parn,'tdiss').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tdiss
          elseif(index(parn,'tice').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tfrmol
          elseif(index(parn,'tfrmol').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tfrmol
          elseif(index(parn,'xc ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xc
          elseif(index(parn,'xsi ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xsi
          elseif(index(parn,'ac ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) ac
          elseif(index(parn,'asi ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) asi
          elseif(index(parn,'kit ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) kit
          elseif(index(parn,'lit ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) lit
          elseif(index(parn,'mit ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) mit
          elseif(index(parn,'ng ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) ng
            ngc = ng/2
          elseif(index(parn,'nn9 ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) nn9
          elseif(index(parn,'nr ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) nr
          elseif(index(parn,'nnt ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) nnt
          elseif(index(parn,'nt ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) nt
          elseif(index(parn,'nrot').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) nrot
          elseif(index(parn,'lrot').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) lrot
          elseif(index(parn,'npl ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) npl
          elseif(index(parn,'nx ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) nx
          elseif(index(parn,'lx ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) lx
          elseif(index(parn,'ngc ').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) ngc
          elseif(index(parn,'moltype').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) moltype
          else
            print '("I don''t understand",a60)', line(1:60)
          endif
        endif
        goto 100

   90   close(unit=12)
        return
      end




      function badnum(val)
!	test val for NaN or Inf

        implicit real*8 (a-h,o-z)
        logical badnum

        if(val.eq.(0.0d+00)) then
          badnum = .false.
        else
          test = val/val
          badnum = (test.ne.test)
        endif

        return
      end




